GLMM_Create_Datasets
2024-04-05
Load packages
# Check you have them and load them
list.of.packages <- c("kableExtra", "tidyr", "ggplot2", "gridExtra", "dplyr", "Hmsc", "jtools", "lubridate", "corrplot", "MuMIn","stringr","sf","raster","leaflet", "tidyverse","htmlwidgets","webshot", "purrr", "magick","forcats","multcomp", "reshape2", "lme4", "tidyr", "stats","RColorBrewer")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
lapply(list.of.packages, require, character.only = TRUE)## Loading required package: kableExtra## Loading required package: tidyr## Loading required package: ggplot2## Loading required package: gridExtra## Loading required package: dplyr##
## Attaching package: 'dplyr'## The following object is masked from 'package:gridExtra':
##
## combine## The following object is masked from 'package:kableExtra':
##
## group_rows## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union## Loading required package: Hmsc## Loading required package: coda## The legacy packages maptools, rgdal, and rgeos, underpinning this package
## will retire shortly. Please refer to R-spatial evolution reports on
## https://r-spatial.org/r/2023/05/15/evolution4.html for details.
## This package is now running under evolution status 0## Loading required package: jtools## Loading required package: lubridate##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union## Loading required package: corrplot## corrplot 0.92 loaded## Loading required package: MuMIn## Loading required package: stringr## Loading required package: sf## Warning: package 'sf' was built under R version 4.3.2## Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUE## Loading required package: raster## Loading required package: sp##
## Attaching package: 'raster'## The following object is masked from 'package:dplyr':
##
## select## Loading required package: leaflet## Loading required package: tidyverse## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ readr 2.1.4
## ✔ purrr 1.0.1 ✔ tibble 3.2.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::combine() masks gridExtra::combine()
## ✖ raster::extract() masks tidyr::extract()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::group_rows() masks kableExtra::group_rows()
## ✖ dplyr::lag() masks stats::lag()
## ✖ raster::select() masks dplyr::select()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## Loading required package: htmlwidgets
##
## Loading required package: webshot
##
## Loading required package: magick
##
## Linking to ImageMagick 6.9.12.98
## Enabled features: cairo, freetype, fftw, ghostscript, heic, lcms, pango, raw, rsvg, webp
## Disabled features: fontconfig, x11
##
## Loading required package: multcomp## Warning: package 'multcomp' was built under R version 4.3.3## Loading required package: mvtnorm## Warning: package 'mvtnorm' was built under R version 4.3.3##
## Attaching package: 'mvtnorm'
##
## The following object is masked from 'package:jtools':
##
## standardize
##
## Loading required package: survival
## Loading required package: TH.data## Warning: package 'TH.data' was built under R version 4.3.3## Loading required package: MASS
##
## Attaching package: 'MASS'
##
## The following objects are masked from 'package:raster':
##
## area, select
##
## The following object is masked from 'package:dplyr':
##
## select
##
##
## Attaching package: 'TH.data'
##
## The following object is masked from 'package:MASS':
##
## geyser
##
## Loading required package: reshape2
##
## Attaching package: 'reshape2'
##
## The following object is masked from 'package:tidyr':
##
## smiths
##
## Loading required package: lme4## Warning: package 'lme4' was built under R version 4.3.2## Loading required package: Matrix## Warning: package 'Matrix' was built under R version 4.3.2##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
##
## Attaching package: 'lme4'
##
## The following object is masked from 'package:raster':
##
## getData
##
## Loading required package: RColorBrewer## [[1]]
## [1] TRUE
##
## [[2]]
## [1] TRUE
##
## [[3]]
## [1] TRUE
##
## [[4]]
## [1] TRUE
##
## [[5]]
## [1] TRUE
##
## [[6]]
## [1] TRUE
##
## [[7]]
## [1] TRUE
##
## [[8]]
## [1] TRUE
##
## [[9]]
## [1] TRUE
##
## [[10]]
## [1] TRUE
##
## [[11]]
## [1] TRUE
##
## [[12]]
## [1] TRUE
##
## [[13]]
## [1] TRUE
##
## [[14]]
## [1] TRUE
##
## [[15]]
## [1] TRUE
##
## [[16]]
## [1] TRUE
##
## [[17]]
## [1] TRUE
##
## [[18]]
## [1] TRUE
##
## [[19]]
## [1] TRUE
##
## [[20]]
## [1] TRUE
##
## [[21]]
## [1] TRUE
##
## [[22]]
## [1] TRUE
##
## [[23]]
## [1] TRUE
##
## [[24]]
## [1] TRUE
##
## [[25]]
## [1] TRUE
##
## [[26]]
## [1] TRUEPrepare Data
locs <- read.csv("Processed_Data/TDN_camera_locations_and_covariates_cleaned.csv", header=T)
# Convert to categorical factors
locs <- locs %>%
mutate_if(is.character,as.factor)
# You should also standardize your covariates - it helps models coverage and facillitates comparison of effects sizes
#NEED TO RUN THIS SO LAT/LONG DONT GET STANDARDIZED
#non response variable species should be scaled but that can manually be done in the GLMM formulas
library(MuMIn)
z_locs <- stdize(locs, omit.cols = c("latitude","longitude","latitude_site","longitude_site","spatial_scale"))Delineate Caribou Seasons
Rather than defining seasons in a community aspect using snowcover (which would be more suited to a multispecies analysis) I want to divide data by biologically relevant “caribou seasons”. The problem here is that these seasons vary by date from year to year and from herd to herd based on GPS collar data. I have 2021-2022 dates for the Ahiak herd from GN reports but I don’t have similar recent seasonal date ranges for the Bathurst/Beverly herds (date ranges used in many reports for the Bathurst herd are from a 2011 study). Do have telemetry data for the Bev/Bath caribou but it’s clipped to the spatial extent of TDN.
Determine first obs, last obs, and number of days observed in TDN for each individual collared caribou using 2021/2022 collar data. Extracted attribute table as csv from shapefile using arcgis
#read in collar data
TDN_Bev_Bath_GPS<-read.csv("Raw_Data/TDN_Bev_Bath_Telemetry.csv")
# Convert LoctnDt to Date format
TDN_Bev_Bath_GPS$LoctnDt <- as.Date(TDN_Bev_Bath_GPS$LoctnDt)Summarize the collar data
# Group by AnimalId, Sex, and Herd, then summarize to find first date, last date, count of observations, and number of unique dates
caribou_indv_summary <- TDN_Bev_Bath_GPS %>%
group_by(AnimlId, SEX, Herd) %>%
summarize(
FirstDate = min(LoctnDt),
LastDate = max(LoctnDt),
NumObservations = n(),
NumUniqueDates = n_distinct(LoctnDt)
)## `summarise()` has grouped output by 'AnimlId', 'SEX'. You can override using
## the `.groups` argument.# Group by Herd, Sex, and then summarize to find first date, last date, count of observations, number of individuals, and number of unique dates
herd_summary <- TDN_Bev_Bath_GPS %>%
group_by(Herd, SEX) %>%
summarize(
FirstDate = min(LoctnDt),
LastDate = max(LoctnDt),
NumObservations = n(),
NumIndividuals = n_distinct(AnimlId),
NumUniqueDates = n_distinct(LoctnDt)
)## `summarise()` has grouped output by 'Herd'. You can override using the
## `.groups` argument.print(herd_summary)## # A tibble: 5 × 7
## # Groups: Herd [3]
## Herd SEX FirstDate LastDate NumObservations NumIndividuals
## <chr> <chr> <date> <date> <int> <int>
## 1 Bathurst Female 2021-10-27 2021-11-15 45 2
## 2 Bathurst Male 2021-10-29 2022-05-06 186 3
## 3 Beverly Female 2021-10-27 2022-01-17 707 16
## 4 Beverly Male 2021-10-29 2022-08-30 1049 12
## 5 Unassigned Female 2021-11-18 2021-11-21 7 1
## # ℹ 1 more variable: NumUniqueDates <int># If you want separate summaries for each Herd irrespective of sex
herd_summary_no_sex <- TDN_Bev_Bath_GPS %>%
group_by(Herd) %>%
summarize(
FirstDate = min(LoctnDt),
LastDate = max(LoctnDt),
NumObservations = n(),
NumIndividuals = n_distinct(AnimlId),
NumUniqueDates = n_distinct(LoctnDt)
)Herd SEX FirstDate LastDate NumObservations NumIndividuals
NumUniqueDates
Looks like the cows are only using TDN during Rut/Breeding, Fall Migration - post-breeding, and the start of Winter according to seasonal cutoff dates other Bathurst papers have used.
Ahiak have slightly different seasonal cutoff dates but collared cows roughly use TDN in 2021/2022 in Late Summer, Fall migration - pre-breeding, Rut/breeding, Fall migration - post-breeding, Winter, no spring migration, then back for Summer and Late Summer.
If we used these dates, how do we account for the different herds (Ahiak winters here but they’re expected to winter on the tundra as opposed to Bev/Bath) and the seasons where we’re less likely to have female caribou (i.e. spring)
Multiseason, DAILY
The telemetry data was informative but still found it difficult to set seasonal cutoff dates based on that alone. Plotting daily detection rate for caribou might be a better way to find natural seasonal breaks and is the most detailed data we have on hand.
Let’s see how this breaks down into 1 day events to help determine seasonal breaks
#Load biweekly detection data
daily_obs<-read.csv("Processed_Data/wildRtrax/TDN_cam_summaries_day.csv")Join this dataframe with the standardized location data:
mod_dat_daily <- left_join(daily_obs, z_locs)## Joining with `by = join_by(location)`And let’s do another raw data check:
boxplot(mod_dat_daily$Caribou~mod_dat_daily$dom_habitat_lcc_hexagon,
las=2,
xlab="lcc landcover",
ylab="Habitat use")
##Add snow based seasons
#format day as Date
mod_dat_daily$day<-as.Date(mod_dat_daily$day)
#create column season
mod_dat_daily <- mod_dat_daily %>%
mutate(snow_season = ifelse(month(day) %in% c(5, 6, 7, 8, 9, 10), "summer", "winter"))
# Get the column names
column_names <- names(mod_dat_daily)
# Reorder the columns
new_column_order <- c(column_names[1:4], "snow_season", column_names[5:length(column_names)])
# Create the reordered data frame
mod_dat_daily <- mod_dat_daily[, new_column_order]
# make it a factor factors
mod_dat_daily <- mod_dat_daily %>%
mutate_if(is.character,as.factor)let’s plot caribou detections summed across locations per day coloured by snow based seasons
# Group by day and summarize the total Caribou count and effort
daily_sum <- mod_dat_daily %>%
group_by(year = year(day), day, snow_season) %>%
summarise(total_caribou = sum(Caribou), total_effort = sum(n_days_effort))## `summarise()` has grouped output by 'year', 'day'. You can override using the
## `.groups` argument.# Plot the data for total Caribou count
ggplot(daily_sum, aes(x = day, y = total_caribou, fill = snow_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Caribou Count Summed Across Locations by Day",
x = "Day",
y = "Total Caribou Count") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8))+
scale_x_date(date_breaks = "1 month", date_labels = "%b %Y")
# Plot the data for capture rate
daily_sum$total_daily_cr <- (daily_sum$total_caribou / daily_sum$total_effort)
ggplot(daily_sum, aes(x = day, y = total_daily_cr, fill = snow_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Caribou Daily CR Across Locations by Day",
x = "Day",
y = "Total Caribou CR") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) +
scale_x_date(date_breaks = "1 month", date_labels = "%b %Y")
Based on the daily detections it looks like caribou detection rate could
be broken down into the following biologically relevant seasons. These
dates roughly conform to date ranges for the Beverly, Bathurst, and
Ahiak herds listed in various GN and GNWT reports.
Typical caribou seasons (Nagy 2011 and GN use 9 seasons): July 29, 2021 - Oct 20, 2021: Summer/Late Summer Oct 21, 2021 - Oct 31, 2021: Pre-breeding Migration Nov 1, 2021 - Dec 14, 2021: Rut/Post-breeding Migration Dec 15, 2021 - Apr 4, 2022: Winter Apr 5, 2022 - May 14, 2022: Spring Migration May 15, 2022 - early June, 2022: Calving early June, 2022 - July 30, 2022: Post-calving/Summer July 31, 2022 - Aug 24, 2022: Late Summer
Simplified 4 seasons (Bathurst Range Plan also uses 4): July 29, 2021 - Oct 20, 2021: Summer Oct 21, 2021 - Dec 14, 2021: Fall Dec 15, 2021 - Apr 4, 2022: Winter Apr 5, 2022 - June 21, 2022: Spring June 22, 2022 - Aug 24, 2022: Summer
##Add simplified 4 seasons
# Define function to assign seasons based on date
assign_season <- function(date) {
case_when(
date >= as.Date("2021-07-29") & date <= as.Date("2021-10-20") ~ "Summer",
date >= as.Date("2021-10-21") & date <= as.Date("2021-12-14") ~ "Fall",
date >= as.Date("2021-12-15") & date <= as.Date("2022-04-04") ~ "Winter",
date >= as.Date("2022-04-05") & date <= as.Date("2022-06-21") ~ "Spring",
date >= as.Date("2022-06-22") & date <= as.Date("2022-08-24") ~ "Summer",
TRUE ~ "Unknown" # Default case
)
}
# Create new column "four_season"
mod_dat_daily <- mod_dat_daily %>%
mutate(four_season = assign_season(day))
# Get the column names
column_names <- names(mod_dat_daily)
# Reorder the columns
new_column_order <- c(column_names[1:4], "snow_season", "four_season", column_names[5:length(column_names)])
# Create the reordered data frame
mod_dat_daily <- mod_dat_daily[, new_column_order]
# Make sure it's a factor
mod_dat_daily <- mod_dat_daily %>%
mutate_if(is.character,as.factor)let’s plot it out
# Group by day and summarize the total Caribou count and effort
daily_sum <- mod_dat_daily %>%
group_by(year = year(day), day, four_season) %>%
summarise(total_caribou = sum(Caribou), total_effort = sum(n_days_effort))## `summarise()` has grouped output by 'year', 'day'. You can override using the
## `.groups` argument.# Plot the data for total Caribou count
ggplot(daily_sum, aes(x = day, y = total_caribou, fill = four_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Caribou Count Summed Across Locations by Day",
x = "Day",
y = "Total Caribou Count") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) +
scale_x_date(date_breaks = "1 month", date_labels = "%b %Y")
# Plot the data for capture rate
daily_sum$total_daily_cr <- (daily_sum$total_caribou / daily_sum$total_effort)
ggplot(daily_sum, aes(x = day, y = total_daily_cr, fill = four_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Caribou Daily CR Across Locations by Day",
x = "Day",
y = "Total Caribou CR") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) +
scale_x_date(date_breaks = "1 month", date_labels = "%b %Y")
###Updated colours for thesis
# Group by day and summarize the total Caribou count and effort
daily_sum <- mod_dat_daily %>%
group_by(year = year(day), day, four_season) %>%
summarise(total_caribou = sum(Caribou), total_effort = sum(n_days_effort))## `summarise()` has grouped output by 'year', 'day'. You can override using the
## `.groups` argument.# Define the colors for each season using specific colors from Paired
season_colors <- c("Winter" = "#A6CEE3", # Light Blue
"Summer" = "#B2DF8A", # Light Green
"Spring" = "#CAB2D6", # Light Purple
"Fall" = "#FDBF6F") # Light Orange
# Plot the data for total Caribou count
ggplot(daily_sum, aes(x = day, y = total_caribou, fill = four_season)) +
geom_bar(stat = "identity", position = "dodge") +
labs(fill = "Season",
x = "Day",
y = "Total Caribou Count") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) +
scale_x_date(date_breaks = "1 month", date_labels = "%b %Y") +
scale_fill_manual(values = season_colors)
# Calculate the total daily capture rate
daily_sum$total_daily_cr <- daily_sum$total_caribou / daily_sum$total_effort
# Plot the data for capture rate
ggplot(daily_sum, aes(x = day, y = total_daily_cr, fill = four_season)) +
geom_bar(stat = "identity", position = "dodge") +
labs(fill = "Season",
x = "Day",
y = "Caribou detection rate") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) +
scale_x_date(date_breaks = "1 month", date_labels = "%b %Y") +
scale_fill_manual(values = season_colors)
##Summary Stats
###Capture Rate by Season
# Summarize the total capture rate by season
seasonal_capture_rate <- daily_sum %>%
group_by(four_season) %>%
summarise(total_daily_cr = sum(total_daily_cr))
print(seasonal_capture_rate)## # A tibble: 4 × 2
## four_season total_daily_cr
## <fct> <dbl>
## 1 Fall 7.28
## 2 Spring 0.421
## 3 Summer 1.32
## 4 Winter 1.77###Effort by season
# Summarize the total effort by season
seasonal_effort <- daily_sum %>%
group_by(four_season) %>%
summarise(total_effort = sum(total_effort))
print(seasonal_effort)## # A tibble: 4 × 2
## four_season total_effort
## <fct> <int>
## 1 Fall 12514
## 2 Spring 21861
## 3 Summer 30992
## 4 Winter 18460#calculate mean daily effort by season
seasonal_mean_effort <- mod_dat_daily %>%
group_by(four_season) %>%
summarise(effort_Mean = mean(n_days_effort))
print(seasonal_mean_effort)## # A tibble: 4 × 2
## four_season effort_Mean
## <fct> <dbl>
## 1 Fall 1
## 2 Spring 1
## 3 Summer 1
## 4 Winter 1###Plot daily effort
# Summarize total days effort for each week
total_days_effort <- mod_dat_daily %>%
group_by(day, four_season) %>%
summarise(total_effort = sum(n_days_effort), .groups = 'drop')
print(total_days_effort)## # A tibble: 392 × 3
## day four_season total_effort
## <date> <fct> <int>
## 1 2021-07-29 Summer 6
## 2 2021-07-30 Summer 6
## 3 2021-07-31 Summer 6
## 4 2021-08-01 Summer 6
## 5 2021-08-02 Summer 6
## 6 2021-08-03 Summer 6
## 7 2021-08-04 Summer 7
## 8 2021-08-05 Summer 7
## 9 2021-08-06 Summer 7
## 10 2021-08-07 Summer 7
## # ℹ 382 more rows# Create the bar plot with colors based on the 'four_season' column
ggplot(total_days_effort, aes(x = day, y = total_effort, fill = four_season)) +
geom_bar(stat = "identity") + # Bar plot with black outlines
scale_fill_manual(values = season_colors) + # Apply the season colors
labs(fill = "Season",
x = "Day of Year",
y = "Total Effort (Camera Days)") +
theme_minimal() +
scale_x_date(date_breaks = "1 month", date_labels = "%b %Y") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
#try out a line plot
ggplot(total_days_effort, aes(x = day, y = total_effort)) +
geom_line(aes(color=four_season)) + # Line plot
geom_point(aes(shape=four_season, color=four_season), size=2) + # Optional: Add points to the line plot
scale_color_manual(values = season_colors) +
labs(color = "Season",
shape = "Season",
x = "Day",
y = "Total Effort (Days)") +
theme_minimal() +
scale_x_date(date_breaks = "1 month", date_labels = "%b %Y") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
Multiseason, WEEKLY detections
# Import the total observations dataset
weekly_obs <- read.csv("Processed_Data/wildRtrax/TDN_cam_summaries_week.csv", header=T)Join this dataframe with the location data, as before:
mod_dat_wk <- left_join(weekly_obs, z_locs)## Joining with `by = join_by(location)`# Remove rows with year 2021 and weeks 30, 31, 32, or 33 due to low survey effort
mod_dat_wk <- mod_dat_wk %>%
filter(!(year == 2021 & week %in% c(30, 31, 32, 33)))Raw data check:
boxplot(mod_dat_wk$Caribou~mod_dat_wk$dom_habitat_lcc_hexagon,
las=2,
xlab="lcc landcover",
ylab="Habitat use")
Add snow based seasons
# Lets create a new column and give it the value summer
mod_dat_wk$snow_season <- "summer"
mod_dat_wk$snow_season[mod_dat_wk$week %in% c(44:53,1:17)] <- "winter"
# Get the column names
column_names <- names(mod_dat_wk)
# Reorder the columns
new_column_order <- c(column_names[1:3], "snow_season", column_names[4:length(column_names)])
# Create the reordered data frame
mod_dat_wk <- mod_dat_wk[, new_column_order]
# make it a factor factors
mod_dat_wk <- mod_dat_wk %>%
mutate_if(is.character,as.factor)Plot caribou detections and rate summed across locations per week
# Convert 'week' and 'year' to a Date object
mod_dat_wk$week_year <- as.Date(paste(mod_dat_wk$year, mod_dat_wk$week, "1", sep="-"), format="%Y-%U-%u")## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid
## Warning in strptime(x, format, tz = "GMT"): (0-based) yday 367 in year 2021 is
## invalid# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk %>%
group_by(week_year, snow_season) %>%
summarise(total_caribou = sum(Caribou),total_effort = sum(n_days_effort))## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.#create capture rate
weekly_sum$total_weekly_cr<-((weekly_sum$total_caribou)/(weekly_sum$total_effort)*7)
# Plot the count data with color based on the 'season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = snow_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Caribou Count Summed Across Region by Week",
x = "Week of the Year",
y = "Total Caribou Count") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8))## Warning: Removed 1 rows containing missing values (`geom_bar()`).
#plot the Capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = snow_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Caribou Weekly CR Across Region by Week",
x = "Week of the Year",
y = "Total Caribou CR") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y.%m.%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8)) ## Warning: Removed 1 rows containing missing values (`geom_bar()`).
let’s repeat this process with the new 4-season caribou based dates from the daily detection data
Should the 4-seasons be based now on week instead of the daily detections even though it’s less detailed?
based on daily detection rate: July 29, 2021 - Oct 20, 2021: Summer Oct 21, 2021 - Dec 14, 2021: Fall Dec 15, 2021 - Apr 4, 2022: Winter Apr 5, 2022 - June 21, 2022: Spring June 22, 2022 - Aug 24, 2022: Summer
based on weekly detection rates? July 26, 2021 - Oct 18, 2021: Summer Oct 25, 2021 - Dec 13, 2021: Fall Dec 20, 2021 - Apr 4, 2022: Winter Apr 11, 2022 - June 21, 2022: Spring June 22, 2022 - Aug 24, 2022: Summer
Add Caribou-centric Seasons
# Let's create a new column for the four seasons based on the specified date ranges
mod_dat_wk$four_season <- "Summer"
mod_dat_wk$four_season[mod_dat_wk$week_year >= as.Date("2021-07-29") & mod_dat_wk$week_year <= as.Date("2021-10-20")] <- "Summer"
mod_dat_wk$four_season[mod_dat_wk$week_year >= as.Date("2021-10-21") & mod_dat_wk$week_year <= as.Date("2021-12-14")] <- "Fall"
mod_dat_wk$four_season[mod_dat_wk$week_year >= as.Date("2021-12-15") & mod_dat_wk$week_year <= as.Date("2022-04-04")] <- "Winter"
mod_dat_wk$four_season[mod_dat_wk$week_year >= as.Date("2022-04-05") & mod_dat_wk$week_year <= as.Date("2022-06-21")] <- "Spring"
mod_dat_wk$four_season[mod_dat_wk$week_year >= as.Date("2022-06-22") & mod_dat_wk$week_year <= as.Date("2022-08-24")] <- "Summer"
# Get the column names
column_names <- names(mod_dat_wk)
# Reorder the columns
new_column_order <- c(column_names[1:3], "four_season", column_names[4:length(column_names)])
# Create the reordered data frame
mod_dat_wk <- mod_dat_wk[, new_column_order]
# Make it a factor
mod_dat_wk <- mod_dat_wk %>%
mutate_if(is.character, as.factor)Plot it out
# Let's look at caribou detections summed across locations per week
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk %>%
group_by(week_year, four_season) %>%
summarise(total_caribou = sum(Caribou), total_effort = sum(n_days_effort))## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.# Create capture rate per 7 days
weekly_sum$total_weekly_cr <- (weekly_sum$total_caribou / weekly_sum$total_effort) * 7
# Plot the count data with color based on the 'four_season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = four_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Caribou Count Summed Across Region by Week",
x = "Week of the Year",
y = "Total Caribou Count") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8))## Warning: Removed 1 rows containing missing values (`geom_bar()`).
# Plot the capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = four_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Caribou Weekly CR Across Region by Week",
x = "Week of the Year",
y = "Total Caribou CR") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) ## Warning: Removed 1 rows containing missing values (`geom_bar()`).
##Summary Stats
###Capture Rate by Season
# Summarize the total capture rate by season
seasonal_capture_rate <- weekly_sum %>%
group_by(four_season) %>%
summarise(total_weekly_cr = sum(total_weekly_cr))
print(seasonal_capture_rate)## # A tibble: 4 × 2
## four_season total_weekly_cr
## <fct> <dbl>
## 1 Fall 7.31
## 2 Spring 0.428
## 3 Summer 1.66
## 4 Winter 1.69# Calculate the average weekly capture rate by season
avg_weekly_cr_season <- weekly_sum %>%
group_by(four_season) %>%
summarise(avg_weekly_cr = mean(total_weekly_cr))
print(avg_weekly_cr_season)## # A tibble: 4 × 2
## four_season avg_weekly_cr
## <fct> <dbl>
## 1 Fall 0.914
## 2 Spring 0.0389
## 3 Summer 0.0875
## 4 Winter 0.106###Effort by season
# Summarize the total effort by season
seasonal_effort_wk <- weekly_sum %>%
group_by(four_season) %>%
summarise(total_effort = sum(total_effort))
print(seasonal_effort_wk)## # A tibble: 4 × 2
## four_season total_effort
## <fct> <int>
## 1 Fall 12490
## 2 Spring 21579
## 3 Summer 30002
## 4 Winter 19238#calculate mean weekly effort by season
seasonal_mean_effort_wk <- mod_dat_wk %>%
group_by(four_season) %>%
summarise(effort_Mean = mean(n_days_effort))
print(seasonal_mean_effort_wk)## # A tibble: 4 × 2
## four_season effort_Mean
## <fct> <dbl>
## 1 Fall 6.21
## 2 Spring 6.84
## 3 Summer 6.49
## 4 Winter 6.21###Effort by Month
# Let's create a new column for the month/year based on the specified date ranges
mod_dat_wk$month[mod_dat_wk$week_year >= as.Date("2021-07-29") & mod_dat_wk$week_year <= as.Date("2021-08-31")] <- "2021-08"
mod_dat_wk$month[mod_dat_wk$week_year >= as.Date("2021-09-01") & mod_dat_wk$week_year <= as.Date("2021-09-30")] <- "2021-09"
mod_dat_wk$month[mod_dat_wk$week_year >= as.Date("2021-10-01") & mod_dat_wk$week_year <= as.Date("2021-10-31")] <- "2021-10"
mod_dat_wk$month[mod_dat_wk$week_year >= as.Date("2021-11-01") & mod_dat_wk$week_year <= as.Date("2021-11-30")] <- "2021-11"
mod_dat_wk$month[mod_dat_wk$week_year >= as.Date("2021-12-01") & mod_dat_wk$week_year <= as.Date("2021-12-31")] <- "2021-12"
mod_dat_wk$month[mod_dat_wk$week_year >= as.Date("2022-01-01") & mod_dat_wk$week_year <= as.Date("2022-01-31")] <- "2022-01"
mod_dat_wk$month[mod_dat_wk$week_year >= as.Date("2022-02-01") & mod_dat_wk$week_year <= as.Date("2022-02-28")] <- "2022-02"
mod_dat_wk$month[mod_dat_wk$week_year >= as.Date("2022-03-01") & mod_dat_wk$week_year <= as.Date("2022-03-31")] <- "2022-03"
mod_dat_wk$month[mod_dat_wk$week_year >= as.Date("2022-04-01") & mod_dat_wk$week_year <= as.Date("2022-04-30")] <- "2022-04"
mod_dat_wk$month[mod_dat_wk$week_year >= as.Date("2022-05-01") & mod_dat_wk$week_year <= as.Date("2022-05-31")] <- "2022-05"
mod_dat_wk$month[mod_dat_wk$week_year >= as.Date("2022-06-01") & mod_dat_wk$week_year <= as.Date("2022-06-30")] <- "2022-06"
mod_dat_wk$month[mod_dat_wk$week_year >= as.Date("2022-07-01") & mod_dat_wk$week_year <= as.Date("2022-07-31")] <- "2022-07"
mod_dat_wk$month[mod_dat_wk$week_year >= as.Date("2022-08-01") & mod_dat_wk$week_year <= as.Date("2022-08-31")] <- "2022-08"
# Get the column names
column_names <- names(mod_dat_wk)
# Reorder the columns
new_column_order <- c(column_names[1:5], "month", column_names[6:length(column_names)])
# Create the reordered data frame
mod_dat_wk <- mod_dat_wk[, new_column_order]
# Make it a factor
mod_dat_wk <- mod_dat_wk %>%
mutate_if(is.character, as.factor)
#calculate mean weekly effort by month/year
monthly_mean_effort_wk <- mod_dat_wk %>%
group_by(month) %>%
summarise(effort_Mean = mean(n_days_effort))
print(monthly_mean_effort_wk)## # A tibble: 14 × 2
## month effort_Mean
## <fct> <dbl>
## 1 2021-08 6.20
## 2 2021-09 6.42
## 3 2021-10 6.83
## 4 2021-11 6.06
## 5 2021-12 5.93
## 6 2022-01 5.57
## 7 2022-02 6.20
## 8 2022-03 6.57
## 9 2022-04 6.87
## 10 2022-05 6.74
## 11 2022-06 7.00
## 12 2022-07 7.00
## 13 2022-08 6.16
## 14 <NA> 1###Capture Rate by Month
# Calculate the average weekly capture rate by month
avg_weekly_cr_month <- mod_dat_wk %>%
group_by(month) %>%
summarise(total_caribou = sum(Caribou), total_effort = sum(n_days_effort)) %>%
mutate(avg_weekly_cr = (total_caribou / total_effort) * 7)
print(avg_weekly_cr_month)## # A tibble: 14 × 4
## month total_caribou total_effort avg_weekly_cr
## <fct> <int> <int> <dbl>
## 1 2021-08 3 1792 0.0117
## 2 2021-09 21 6074 0.0242
## 3 2021-10 299 8124 0.258
## 4 2021-11 1164 8058 1.01
## 5 2021-12 297 4186 0.497
## 6 2022-01 53 3742 0.0991
## 7 2022-02 56 4412 0.0888
## 8 2022-03 14 7252 0.0135
## 9 2022-04 60 7853 0.0535
## 10 2022-05 56 9647 0.0406
## 11 2022-06 11 8094 0.00951
## 12 2022-07 16 8034 0.0139
## 13 2022-08 294 5946 0.346
## 14 <NA> 4 95 0.295###Plot weekly effort
# Summarize total days effort for each week
total_weeks_effort <- mod_dat_wk %>%
group_by(week_year, four_season) %>%
summarise(total_effort = sum(n_days_effort), .groups = 'drop')
# Convert 'week_year' to Date format if not already done
total_weeks_effort$week_year <- as.Date(total_weeks_effort$week_year, format = "%Y-%m-%d")
# Create the bar plot with colors based on the 'four_season' column
ggplot(total_weeks_effort, aes(x = week_year, y = total_effort, fill = four_season)) +
geom_bar(stat = "identity") + # Bar plot with black outlines
scale_fill_manual(values = season_colors) + # Apply the season colors
labs(fill = "Season",
x = "Week of Year",
y = "Total Effort (Camera Days)") +
theme_minimal() +
scale_x_date(date_breaks = "1 month", date_labels = "%b %Y") + # Format the x-axis for weekly data
theme(axis.text.x = element_text(angle = 45, hjust = 1))## Warning: Removed 1 rows containing missing values (`position_stack()`).
##Consolidate Ungulates and Predators
# Create new columns "Ungulates" and "Predators" by summing the respective columns
mod_dat_wk$Ungulates <- mod_dat_wk$Muskox + mod_dat_wk$Moose
mod_dat_wk$Predators <- mod_dat_wk$Gray.Wolf + mod_dat_wk$Grizzly.Bear##Plotting Create the same plots for each of the other focal species Muskox
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk %>%
group_by(week_year, snow_season) %>%
summarise(total_caribou = sum(Muskox),total_effort = sum(n_days_effort))## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.#create capture rate
weekly_sum$total_weekly_cr<-((weekly_sum$total_caribou)/(weekly_sum$total_effort)*7)
# Plot the count data with color based on the 'season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = snow_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Muskox Count Summed Across Region by Week",
x = "Week of the Year",
y = "Total Muskox Count") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8))## Warning: Removed 1 rows containing missing values (`geom_bar()`).
#plot the Capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = snow_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Muskox Weekly CR Across Region by Week",
x = "Week of the Year",
y = "Total Muskox CR") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y.%m.%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8)) ## Warning: Removed 1 rows containing missing values (`geom_bar()`).
# Let's look at caribou detections summed across locations per week
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk %>%
group_by(week_year, four_season) %>%
summarise(total_caribou = sum(Muskox), total_effort = sum(n_days_effort))## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.# Create capture rate
weekly_sum$total_weekly_cr <- (weekly_sum$total_caribou / weekly_sum$total_effort) * 7
# Plot the count data with color based on the 'four_season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = four_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Muskox Count Summed Across Region by Week",
x = "Week of the Year",
y = "Total Muskox Count") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8))## Warning: Removed 1 rows containing missing values (`geom_bar()`).
# Plot the capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = four_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Muskox Weekly CR Across Region by Week",
x = "Week of the Year",
y = "Total Muskox CR") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) ## Warning: Removed 1 rows containing missing values (`geom_bar()`).
Moose
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk %>%
group_by(week_year, snow_season) %>%
summarise(total_caribou = sum(Moose),total_effort = sum(n_days_effort))## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.#create capture rate
weekly_sum$total_weekly_cr<-((weekly_sum$total_caribou)/(weekly_sum$total_effort)*7)
# Plot the count data with color based on the 'season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = snow_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Moose Count Summed Across Region by Week",
x = "Week of the Year",
y = "Total Moose Count") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8))## Warning: Removed 1 rows containing missing values (`geom_bar()`).
#plot the Capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = snow_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Moose Weekly CR Across Region by Week",
x = "Week of the Year",
y = "Total Moose CR") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y.%m.%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8)) ## Warning: Removed 1 rows containing missing values (`geom_bar()`).
# Let's look at caribou detections summed across locations per week
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk %>%
group_by(week_year, four_season) %>%
summarise(total_caribou = sum(Moose), total_effort = sum(n_days_effort))## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.# Create capture rate
weekly_sum$total_weekly_cr <- (weekly_sum$total_caribou / weekly_sum$total_effort) * 7
# Plot the count data with color based on the 'four_season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = four_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Moose Count Summed Across Region by Week",
x = "Week of the Year",
y = "Total Moose Count") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8))## Warning: Removed 1 rows containing missing values (`geom_bar()`).
# Plot the capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = four_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Moose Weekly CR Across Region by Week",
x = "Week of the Year",
y = "Total Moose CR") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) ## Warning: Removed 1 rows containing missing values (`geom_bar()`).
Gray.Wolf
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk %>%
group_by(week_year, snow_season) %>%
summarise(total_caribou = sum(Gray.Wolf),total_effort = sum(n_days_effort))## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.#create capture rate
weekly_sum$total_weekly_cr<-((weekly_sum$total_caribou)/(weekly_sum$total_effort)*7)
# Plot the count data with color based on the 'season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = snow_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Gray Wolf Count Summed Across Region by Week",
x = "Week of the Year",
y = "Total Gray Wolf Count") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8))## Warning: Removed 1 rows containing missing values (`geom_bar()`).
#plot the Capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = snow_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Gray Wolf Weekly CR Across Region by Week",
x = "Week of the Year",
y = "Total Gray Wolf CR") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y.%m.%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8)) ## Warning: Removed 1 rows containing missing values (`geom_bar()`).
# Let's look at caribou detections summed across locations per week
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk %>%
group_by(week_year, four_season) %>%
summarise(total_caribou = sum(Gray.Wolf), total_effort = sum(n_days_effort))## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.# Create capture rate
weekly_sum$total_weekly_cr <- (weekly_sum$total_caribou / weekly_sum$total_effort) * 7
# Plot the count data with color based on the 'four_season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = four_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Gray Wolf Count Summed Across Region by Week",
x = "Week of the Year",
y = "Total Gray Wolf Count") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8))## Warning: Removed 1 rows containing missing values (`geom_bar()`).
# Plot the capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = four_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Gray Wolf Weekly CR Across Region by Week",
x = "Week of the Year",
y = "Total Gray Wolf CR") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) ## Warning: Removed 1 rows containing missing values (`geom_bar()`).
Grizzly.Bear
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk %>%
group_by(week_year, snow_season) %>%
summarise(total_caribou = sum(Grizzly.Bear),total_effort = sum(n_days_effort))## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.#create capture rate
weekly_sum$total_weekly_cr<-((weekly_sum$total_caribou)/(weekly_sum$total_effort)*7)
# Plot the count data with color based on the 'season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = snow_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Grizzly Bear Count Summed Across Region by Week",
x = "Week of the Year",
y = "Total Grizzly Bear Count") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8))## Warning: Removed 1 rows containing missing values (`geom_bar()`).
#plot the Capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = snow_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Grizzly Bear Weekly CR Across Region by Week",
x = "Week of the Year",
y = "Total Grizzly Bear CR") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y.%m.%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8)) ## Warning: Removed 1 rows containing missing values (`geom_bar()`).
# Let's look at caribou detections summed across locations per week
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk %>%
group_by(week_year, four_season) %>%
summarise(total_caribou = sum(Grizzly.Bear), total_effort = sum(n_days_effort))## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.# Create capture rate
weekly_sum$total_weekly_cr <- (weekly_sum$total_caribou / weekly_sum$total_effort) * 7
# Plot the count data with color based on the 'four_season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = four_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Grizzly Bear Count Summed Across Region by Week",
x = "Week of the Year",
y = "Total Grizzly Bear Count") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8))## Warning: Removed 1 rows containing missing values (`geom_bar()`).
# Plot the capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = four_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Grizzly Bear Weekly CR Across Region by Week",
x = "Week of the Year",
y = "Total Grizzly Bear CR") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) ## Warning: Removed 1 rows containing missing values (`geom_bar()`).
Ungulates
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk %>%
group_by(week_year, snow_season) %>%
summarise(total_caribou = sum(Ungulates),total_effort = sum(n_days_effort))## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.#create capture rate
weekly_sum$total_weekly_cr<-((weekly_sum$total_caribou)/(weekly_sum$total_effort)*7)
# Plot the count data with color based on the 'season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = snow_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Ungulate Count Summed Across Region by Week",
x = "Week of the Year",
y = "Total Ungulate Count") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8))## Warning: Removed 1 rows containing missing values (`geom_bar()`).
#plot the Capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = snow_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Ungulate Weekly CR Across Region by Week",
x = "Week of the Year",
y = "Total Ungulate CR") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y.%m.%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8)) ## Warning: Removed 1 rows containing missing values (`geom_bar()`).
# Let's look at caribou detections summed across locations per week
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk %>%
group_by(week_year, four_season) %>%
summarise(total_caribou = sum(Ungulates), total_effort = sum(n_days_effort))## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.# Create capture rate
weekly_sum$total_weekly_cr <- (weekly_sum$total_caribou / weekly_sum$total_effort) * 7
# Plot the count data with color based on the 'four_season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = four_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Ungulate Count Summed Across Region by Week",
x = "Week of the Year",
y = "Total Ungulate Count") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8))## Warning: Removed 1 rows containing missing values (`geom_bar()`).
# Plot the capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = four_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Ungulate Weekly CR Across Region by Week",
x = "Week of the Year",
y = "Total Ungulate CR") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) ## Warning: Removed 1 rows containing missing values (`geom_bar()`).
Predators
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk %>%
group_by(week_year, snow_season) %>%
summarise(total_caribou = sum(Predators),total_effort = sum(n_days_effort))## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.#create capture rate
weekly_sum$total_weekly_cr<-((weekly_sum$total_caribou)/(weekly_sum$total_effort)*7)
# Plot the count data with color based on the 'season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = snow_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Predator Count Summed Across Region by Week",
x = "Week of the Year",
y = "Total Predator Count") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8))## Warning: Removed 1 rows containing missing values (`geom_bar()`).
#plot the Capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = snow_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Predator Weekly CR Across Region by Week",
x = "Week of the Year",
y = "Total Predator CR") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y.%m.%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8)) ## Warning: Removed 1 rows containing missing values (`geom_bar()`).
# Let's look at caribou detections summed across locations per week
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk %>%
group_by(week_year, four_season) %>%
summarise(total_caribou = sum(Predators), total_effort = sum(n_days_effort))## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.# Create capture rate
weekly_sum$total_weekly_cr <- (weekly_sum$total_caribou / weekly_sum$total_effort) * 7
# Plot the count data with color based on the 'four_season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = four_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Predator Count Summed Across Region by Week",
x = "Week of the Year",
y = "Total Predator Count") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8))## Warning: Removed 1 rows containing missing values (`geom_bar()`).
# Plot the capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = four_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Predator Weekly CR Across Region by Week",
x = "Week of the Year",
y = "Total Predator CR") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) ## Warning: Removed 1 rows containing missing values (`geom_bar()`).
map it out
# Aggregate data by location and week to get total caribou detections
#total_caribou_by_location_week <- mod_dat_wk %>%
# group_by(location, week) %>%
#summarise(total_caribou = sum(Caribou))
# Merge total caribou detections with latitude and longitude data
#locations_week <- mod_dat_wk %>%
# distinct(location, latitude, longitude)
#locations_week <- left_join(locations_week, total_caribou_by_location_week, by = "location") %>%
# mutate(has_detections = if_else(is.na(total_caribou) | total_caribou == 0, FALSE, TRUE))
# Create leaflet maps for each week
#maps_folder <- "maps"
#dir.create(maps_folder, showWarnings = FALSE)
#for (week in unique(mod_dat_wk$week)) {
# map_data_week <- filter(locations_week, week == !!week)
#map <- leaflet(map_data_week) %>%
# addTiles() %>%
#addCircleMarkers(
# ~longitude,
#~latitude,
#radius = ~if_else(has_detections, sqrt(total_caribou) * 2, 2),
#color = ~if_else(has_detections, "blue", "grey"),
#popup = ~paste("Location:", location, "<br>", "Total Caribou Detections:", total_caribou)
#)
# Save the map as HTML
#map_file <- file.path(maps_folder, paste0("map_week_", week, ".html"))
# saveWidget(map, map_file, selfcontained = TRUE)
#}# Function to read HTML files and convert them to images
#html_to_image <- function(html_file) {
# webshot::webshot(html_file, file = tempfile(fileext = ".png"))
#}
# Read each HTML file, convert to images, and store them in a list
#map_images <- lapply(list.files("maps", pattern = "*.html", full.names = TRUE), html_to_image)
# Create the animated GIF
#gif_file <- "caribou_detections.gif"
#image_write(image_join(map_images), gif_file)
# Display the GIF
#browseURL(gif_file)##Target Analysis
# Reorder the levels of the target_cons column
mod_dat_wk$target_cons <- fct_relevel(mod_dat_wk$target_cons, "None")
#check the levels were reordered correctly. None should be first to act as the control feature.
levels(mod_dat_wk$target_cons)## [1] "None" "Esker" "Shoreline" "Trail" "Wetland"# Group by target_cons and count unique locations within each group
unique_locations_by_target <- mod_dat_wk %>%
group_by(target_cons) %>%
summarise(unique_locations = n_distinct(location))
# Print the result
print(unique_locations_by_target)## # A tibble: 5 × 2
## target_cons unique_locations
## <fct> <int>
## 1 None 218
## 2 Esker 7
## 3 Shoreline 20
## 4 Trail 27
## 5 Wetland 35#run anova to determine whether targets are actually significantly related to our response variable (Caribou occurrences in a given week at a given camera)
region_wk_aov <- aov(Caribou ~ target_cons, data=mod_dat_wk)
summary(region_wk_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_cons 4 50 12.62 9.418 1.34e-07 ***
## Residuals 12882 17260 1.34
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_cons, data = mod_dat_wk)
##
## $target_cons
## diff lwr upr p adj
## Esker-None 0.33851816 0.07685521 0.600181105 0.0038263
## Shoreline-None 0.05382748 -0.06082899 0.168483953 0.7031152
## Trail-None -0.15620948 -0.25213604 -0.060282928 0.0000876
## Wetland-None -0.04684874 -0.13440328 0.040705800 0.5887736
## Shoreline-Esker -0.28469068 -0.56653803 -0.002843318 0.0463407
## Trail-Esker -0.49472764 -0.76948887 -0.219966416 0.0000090
## Wetland-Esker -0.38536690 -0.65731841 -0.113415387 0.0010530
## Trail-Shoreline -0.21003696 -0.35206633 -0.068007604 0.0005275
## Wetland-Shoreline -0.10067622 -0.23719081 0.035838367 0.2601690
## Wetland-Trail 0.10936074 -0.01184994 0.230571428 0.0995992perform the same tests on the new target_hydro (wetland shoreline) and target_linear (trail esker) categories. These targets have been separated because their predicted to have similar effects on caribou occurrences.
# Reorder the levels of the target_hydro column
mod_dat_wk$target_hydro <- fct_relevel(mod_dat_wk$target_hydro, "None")
#check the levels were reordered correctly. None should be first to act as the control feature.
levels(mod_dat_wk$target_hydro)## [1] "None" "Shoreline" "Wetland"#run anova on hydro targets
region_wk_hydro_aov <- aov(Caribou ~ target_hydro, data=mod_dat_wk)
summary(region_wk_hydro_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_hydro 2 5 2.716 2.022 0.132
## Residuals 12884 17305 1.343# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_hydro_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_hydro, data = mod_dat_wk)
##
## $target_hydro
## diff lwr upr p adj
## Shoreline-None 0.06728355 -0.03081301 0.16538011 0.2423553
## Wetland-None -0.03339267 -0.10800855 0.04122321 0.5459357
## Wetland-Shoreline -0.10067622 -0.21811018 0.01675774 0.1099318Not significant.
# Reorder the levels of the target_linear column
mod_dat_wk$target_linear <- fct_relevel(mod_dat_wk$target_linear, "None")
#check the levels were reordered correctly. None should be first to act as the control feature.
levels(mod_dat_wk$target_linear)## [1] "None" "Esker" "Trail"#run anova on linear targets
region_wk_linear_aov <- aov(Caribou ~ target_linear, data=mod_dat_wk)
summary(region_wk_linear_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_linear 2 45 22.41 16.72 5.6e-08 ***
## Residuals 12884 17265 1.34
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_linear_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_linear, data = mod_dat_wk)
##
## $target_linear
## diff lwr upr p adj
## Esker-None 0.3408211 0.1163509 0.56529136 0.0010882
## Trail-None -0.1539065 -0.2353337 -0.07247934 0.0000282
## Trail-Esker -0.4947276 -0.7308164 -0.25863884 0.0000027 Df Sum Sq Mean Sq F value Pr(>F) target_linear 2 45 22.351 16.79 5.23e-08 ***
$target_linear diff lwr upr p adj Esker-None 0.3412555 0.1175085 0.56500256 0.0010243 Trail-None -0.1534721 -0.2346369 -0.07230728 0.0000280
Summer, WEEKLY detection data
# Create a new data frame for summer
mod_dat_wk_summer <- mod_dat_wk[mod_dat_wk$four_season == "Summer", ]let’s try to map er out
# Aggregate data by location to get total caribou detections
total_caribou_by_location <- mod_dat_wk_summer %>%
group_by(location) %>%
summarise(total_caribou = sum(Caribou))
# Merge total caribou detections with latitude and longitude data
locations <- mod_dat_wk_summer %>%
distinct(location, latitude, longitude)
locations <- left_join(locations, total_caribou_by_location, by = "location") %>%
mutate(has_detections = if_else(is.na(total_caribou) | total_caribou == 0, FALSE, TRUE))
# Create leaflet map
leaflet(locations) %>%
addTiles() %>%
addCircleMarkers(
~longitude,
~latitude,
radius = ~if_else(has_detections, sqrt(total_caribou) * 2, 2),
color = ~if_else(has_detections, "blue", "grey"),
popup = ~paste("Location:", location, "<br>", "Total Caribou Detections:", total_caribou)
)##Target analysis
#run anova to determine whether targets are actually significantly related to our response variable (Caribou occurrences in a given week at a given camera)
region_wk_smr_aov <- aov(Caribou ~ target_cons, data=mod_dat_wk_summer)
summary(region_wk_smr_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_cons 4 10 2.5275 3.583 0.00636 **
## Residuals 4619 3258 0.7053
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_smr_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_cons, data = mod_dat_wk_summer)
##
## $target_cons
## diff lwr upr p adj
## Esker-None 0.06398224 -2.504331e-01 0.37839758 0.9813072
## Shoreline-None 0.13661331 7.223356e-05 0.27315439 0.0498024
## Trail-None -0.05386288 -1.714791e-01 0.06375331 0.7220456
## Wetland-None -0.07473195 -1.819719e-01 0.03250796 0.3164885
## Shoreline-Esker 0.07263107 -2.654815e-01 0.41074361 0.9771589
## Trail-Esker -0.11784512 -4.487684e-01 0.21307812 0.8679521
## Wetland-Esker -0.13871419 -4.660932e-01 0.18866483 0.7762729
## Trail-Shoreline -0.19047619 -3.616388e-01 -0.01931355 0.0203240
## Wetland-Shoreline -0.21134526 -3.755509e-01 -0.04713967 0.0040974
## Wetland-Trail -0.02086907 -1.697099e-01 0.12797172 0.9954618#run anova on hydro targets
region_wk_smr_hydro_aov <- aov(Caribou ~ target_hydro, data=mod_dat_wk_summer)
summary(region_wk_smr_hydro_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_hydro 2 9 4.373 6.2 0.00205 **
## Residuals 4621 3259 0.705
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_smr_hydro_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_hydro, data = mod_dat_wk_summer)
##
## $target_hydro
## diff lwr upr p adj
## Shoreline-None 0.14180405 0.02513743 0.25847067 0.0122279
## Wetland-None -0.06954121 -0.16085804 0.02177562 0.1745725
## Wetland-Shoreline -0.21134526 -0.35241929 -0.07027123 0.0013034 Df Sum Sq Mean Sq F value Pr(>F) target_hydro 2 9 4.376 6.301 0.00185 **
diff lwr upr p adjShoreline-None 0.14233231 0.02655834 0.25810628 0.0110639
#run anova on linear targets
region_wk_smr_linear_aov <- aov(Caribou ~ target_linear, data=mod_dat_wk_summer)
summary(region_wk_smr_linear_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_linear 2 1 0.7016 0.993 0.371
## Residuals 4621 3267 0.7069# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_smr_linear_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_linear, data = mod_dat_wk_summer)
##
## $target_linear
## diff lwr upr p adj
## Esker-None 0.06338601 -0.2066014 0.33337343 0.8462692
## Trail-None -0.05445910 -0.1544370 0.04551876 0.4082263
## Trail-Esker -0.11784512 -0.4024715 0.16678121 0.5954731Not significant
Fall, WEEKLY detection data
# Create a new data frame for winter
mod_dat_wk_fall <- mod_dat_wk[mod_dat_wk$four_season == "Fall", ]plot it out
# Aggregate data by location to get total caribou detections
total_caribou_by_location <- mod_dat_wk_fall %>%
group_by(location) %>%
summarise(total_caribou = sum(Caribou))
# Merge total caribou detections with latitude and longitude data
locations <- mod_dat_wk_fall %>%
distinct(location, latitude, longitude)
locations <- left_join(locations, total_caribou_by_location, by = "location") %>%
mutate(has_detections = if_else(is.na(total_caribou) | total_caribou == 0, FALSE, TRUE))
# Create leaflet map
leaflet(locations) %>%
addTiles() %>%
addCircleMarkers(
~longitude,
~latitude,
radius = ~if_else(has_detections, sqrt(total_caribou) * 2, 2),
color = ~if_else(has_detections, "blue", "grey"),
popup = ~paste("Location:", location, "<br>", "Total Caribou Detections:", total_caribou)
)##Target analysis
#run anova to determine whether targets are actually significantly related to our response variable (Caribou occurrences in a given week at a given camera)
region_wk_fll_aov <- aov(Caribou ~ target_cons, data=mod_dat_wk_fall)
summary(region_wk_fll_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_cons 4 126 31.513 5.409 0.000249 ***
## Residuals 2007 11694 5.826
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_fll_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_cons, data = mod_dat_wk_fall)
##
## $target_cons
## diff lwr upr p adj
## Esker-None 0.9786811 -0.301492928 2.25885517 0.2259125
## Shoreline-None 0.1107002 -0.506163770 0.72756424 0.9883063
## Trail-None -0.7448285 -1.257196295 -0.23246068 0.0007122
## Wetland-None -0.1032963 -0.563718174 0.35712566 0.9731459
## Shoreline-Esker -0.8679809 -2.267564708 0.53160294 0.4383318
## Trail-Esker -1.7235096 -3.080283150 -0.36673606 0.0048540
## Wetland-Esker -1.0819774 -2.419999057 0.25604431 0.1771240
## Trail-Shoreline -0.8555287 -1.618752187 -0.09230526 0.0189988
## Wetland-Shoreline -0.2139965 -0.943364456 0.51537147 0.9303287
## Wetland-Trail 0.6415322 -0.001878368 1.28494283 0.0511013#run anova on hydro targets
region_wk_fll_hydro_aov <- aov(Caribou ~ target_hydro, data=mod_dat_wk_fall)
summary(region_wk_fll_hydro_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_hydro 2 4 2.087 0.355 0.701
## Residuals 2009 11815 5.881# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_fll_hydro_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_hydro, data = mod_dat_wk_fall)
##
## $target_hydro
## diff lwr upr p adj
## Shoreline-None 0.17914083 -0.3505139 0.7087955 0.7072086
## Wetland-None -0.03485566 -0.4285429 0.3588316 0.9765087
## Wetland-Shoreline -0.21399649 -0.8435072 0.4155143 0.7047199Not significant
#run anova on linear targets
region_wk_fll_linear_aov <- aov(Caribou ~ target_linear, data=mod_dat_wk_fall)
summary(region_wk_fll_linear_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_linear 2 122 61.02 10.48 2.96e-05 ***
## Residuals 2009 11698 5.82
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_fll_linear_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_linear, data = mod_dat_wk_fall)
##
## $target_linear
## diff lwr upr p adj
## Esker-None 0.9847773 -0.1125523 2.0821070 0.0890824
## Trail-None -0.7387323 -1.1735898 -0.3038747 0.0002068
## Trail-Esker -1.7235096 -2.8886724 -0.5583468 0.0015442 Df Sum Sq Mean Sq F value Pr(>F) target_linear 2 122 60.97 10.51 2.87e-05 ***
diff lwr upr p adjTrail-None -0.7381761 -1.1721914 -0.3041608 0.0002027
Winter, WEEKLY detection data
# Create a new data frame for winter
mod_dat_wk_winter <- mod_dat_wk[mod_dat_wk$four_season == "Winter", ]plot it out
# Aggregate data by location to get total caribou detections
total_caribou_by_location <- mod_dat_wk_winter %>%
group_by(location) %>%
summarise(total_caribou = sum(Caribou))
# Merge total caribou detections with latitude and longitude data
locations <- mod_dat_wk_winter %>%
distinct(location, latitude, longitude)
locations <- left_join(locations, total_caribou_by_location, by = "location") %>%
mutate(has_detections = if_else(is.na(total_caribou) | total_caribou == 0, FALSE, TRUE))
# Create leaflet map
leaflet(locations) %>%
addTiles() %>%
addCircleMarkers(
~longitude,
~latitude,
radius = ~if_else(has_detections, sqrt(total_caribou) * 2, 2),
color = ~if_else(has_detections, "blue", "grey"),
popup = ~paste("Location:", location, "<br>", "Total Caribou Detections:", total_caribou)
)##Target analysis
#run anova to determine whether targets are actually significantly related to our response variable (Caribou occurrences in a given week at a given camera)
region_wk_wtr_aov <- aov(Caribou ~ target_cons, data=mod_dat_wk_winter)
summary(region_wk_wtr_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_cons 4 9.6 2.3953 12.29 6.41e-10 ***
## Residuals 3091 602.2 0.1948
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_wtr_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_cons, data = mod_dat_wk_winter)
##
## $target_cons
## diff lwr upr p adj
## Esker-None 0.448791822 0.25654490 0.64103875 0.0000000
## Shoreline-None -0.060501372 -0.15145889 0.03045614 0.3646771
## Trail-None -0.051592794 -0.12328805 0.02010247 0.2839604
## Wetland-None -0.004043230 -0.07048943 0.06240297 0.9998297
## Shoreline-Esker -0.509293194 -0.71877657 -0.29980982 0.0000000
## Trail-Esker -0.500384615 -0.70225147 -0.29851776 0.0000000
## Wetland-Esker -0.452835052 -0.65289782 -0.25277228 0.0000000
## Trail-Shoreline 0.008908578 -0.10093047 0.11874762 0.9994679
## Wetland-Shoreline 0.056458142 -0.05002896 0.16294524 0.5971136
## Wetland-Trail 0.047549564 -0.04303986 0.13813898 0.6065337#run anova on hydro targets
region_wk_wtr_hydro_aov <- aov(Caribou ~ target_hydro, data=mod_dat_wk_winter)
summary(region_wk_wtr_hydro_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_hydro 2 0.7 0.3300 1.67 0.188
## Residuals 3093 611.1 0.1976# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_wtr_hydro_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_hydro, data = mod_dat_wk_winter)
##
## $target_hydro
## diff lwr upr p adj
## Shoreline-None -0.060971779 -0.13919861 0.01725505 0.1607324
## Wetland-None -0.004513637 -0.06136034 0.05233307 0.9810722
## Wetland-Shoreline 0.056458142 -0.03567087 0.14858715 0.3220188Not significant
#run anova on linear targets
region_wk_wtr_linear_aov <- aov(Caribou ~ target_linear, data=mod_dat_wk_winter)
summary(region_wk_wtr_linear_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_linear 2 8.9 4.469 22.93 1.3e-10 ***
## Residuals 3093 602.9 0.195
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_wtr_linear_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_linear, data = mod_dat_wk_winter)
##
## $target_linear
## diff lwr upr p adj
## Esker-None 0.45359758 0.2887224 0.61847281 0.0000000
## Trail-None -0.04678703 -0.1075309 0.01395679 0.1677186
## Trail-Esker -0.50038462 -0.6738460 -0.32692325 0.0000000 Df Sum Sq Mean Sq F value Pr(>F) target_linear 2 8.9 4.470 23.12 1.08e-10 ***
diff lwr upr p adjEsker-None 0.45396375 0.2897541 0.61817339 0.0000000
Spring, WEEKLY detection data
# Create a new data frame for winter
mod_dat_wk_spri <- mod_dat_wk[mod_dat_wk$four_season == "Spring", ]plot it out
# Aggregate data by location to get total caribou detections
total_caribou_by_location <- mod_dat_wk_spri %>%
group_by(location) %>%
summarise(total_caribou = sum(Caribou))
# Merge total caribou detections with latitude and longitude data
locations <- mod_dat_wk_spri %>%
distinct(location, latitude, longitude)
locations <- left_join(locations, total_caribou_by_location, by = "location") %>%
mutate(has_detections = if_else(is.na(total_caribou) | total_caribou == 0, FALSE, TRUE))
# Create leaflet map
leaflet(locations) %>%
addTiles() %>%
addCircleMarkers(
~longitude,
~latitude,
radius = ~if_else(has_detections, sqrt(total_caribou) * 2, 2),
color = ~if_else(has_detections, "blue", "grey"),
popup = ~paste("Location:", location, "<br>", "Total Caribou Detections:", total_caribou)
)##Target analysis
#run anova to determine whether targets are actually significantly related to our response variable (Caribou occurrences in a given week at a given camera)
region_wk_spr_aov <- aov(Caribou ~ target_cons, data=mod_dat_wk_spri)
summary(region_wk_spr_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_cons 4 0.5 0.1283 0.61 0.655
## Residuals 3150 662.0 0.2102# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_spr_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_cons, data = mod_dat_wk_spri)
##
## $target_cons
## diff lwr upr p adj
## Esker-None -0.04172156 -0.28394226 0.20049913 0.9900043
## Shoreline-None 0.02200393 -0.06943901 0.11344687 0.9654090
## Trail-None -0.03130490 -0.10955760 0.04694780 0.8108308
## Wetland-None -0.01943744 -0.09048961 0.05161473 0.9454063
## Shoreline-Esker 0.06372549 -0.19251186 0.31996284 0.9610529
## Trail-Esker 0.01041667 -0.24141513 0.26224846 0.9999634
## Wetland-Esker 0.02228412 -0.22740404 0.27197228 0.9992236
## Trail-Shoreline -0.05330882 -0.16780860 0.06119095 0.7092830
## Wetland-Shoreline -0.04144137 -0.15114608 0.06826335 0.8410457
## Wetland-Trail 0.01186746 -0.08711132 0.11084623 0.9975265#run anova on hydro targets
region_wk_spr_hydro_aov <- aov(Caribou ~ target_hydro, data=mod_dat_wk_spri)
summary(region_wk_spr_hydro_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_hydro 2 0.2 0.1118 0.532 0.588
## Residuals 3152 662.3 0.2101# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_spr_hydro_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_hydro, data = mod_dat_wk_spri)
##
## $target_hydro
## diff lwr upr p adj
## Shoreline-None 0.02591685 -0.05224127 0.10407497 0.7168795
## Wetland-None -0.01552452 -0.07605275 0.04500371 0.8193488
## Wetland-Shoreline -0.04144137 -0.13568022 0.05279749 0.5572306Not Significant
#run anova on linear targets
region_wk_spr_linear_aov <- aov(Caribou ~ target_linear, data=mod_dat_wk_spri)
summary(region_wk_spr_linear_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_linear 2 0.3 0.1404 0.668 0.513
## Residuals 3152 662.2 0.2101# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_spr_linear_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_linear, data = mod_dat_wk_spri)
##
## $target_linear
## diff lwr upr p adj
## Esker-None -0.04084507 -0.24866735 0.16697721 0.8895553
## Trail-None -0.03042840 -0.09689405 0.03603725 0.5306441
## Trail-Esker 0.01041667 -0.20590322 0.22673656 0.9929959Not Significant
Caribou Sites Only
Subset detection data to only the clusters which captured caribou at a least one of their cameras. Must be subset after joining with covariates to filter by location_site. Must decide on response variable time (i.e. week, biweek) before trying this to cut down on time.
I think weekly is probably the way to go. it’s a good balance of high detail without being so coarse as to become irrelevant for species co-occurrences. It also provides adequate survey replicas when running 4 separate seasonal models.
##UPDATE: 2024-03-21 ##############################################################
Realized that you can’t just subset the data based on caribou detections being greater than 1 because you want to keep ALL cameras and ALL sampling periods at clusters with at least one caribou, whereas the original code was only taking the cameras/sampling periods with caribou and dumping everything else. Need to figure out how to keep only those cameras (i.e. location) which had at least 1 caribou detection at the cluster level (i.e. location_site).
##UPDATE: 2024-03-22 ##############################################################
OKAY I THINK I FIXED IT.
weekly caribou data
# Filter rows where "Caribou" is greater than or equal to 1
caribou_data <- mod_dat_wk %>%
filter(Caribou >= 1) %>%
dplyr::select(location, location_site) %>%
distinct()
# Filter mod_dat_wk for rows with locations that have caribou detections or share the same location_site
mod_dat_wk_cbou <- mod_dat_wk %>%
filter(location %in% caribou_data$location | location_site %in% caribou_data$location_site)Better to use biological rationale for choosing the local sites rather than a detection based criteria. Proportion of “TUNDRA” vs “TAIGA”? Compare with locations from caribou based data above
# Create a new dataframe with rows where z.TUNDRA > z.TAIGA
mod_dat_wk_tundra <- mod_dat_wk %>%
filter(z.TUNDRA > z.TAIGA)Compare with locations from caribou detection filtered data
# Get unique locations in mod_dat_wk_tundra
unique_locations_tundra <- unique(mod_dat_wk_tundra$location)
# Get unique locations in mod_dat_wk_cbou
unique_locations_cbou <- unique(mod_dat_wk_cbou$location)
# Compare unique locations
common_locations <- intersect(unique_locations_tundra, unique_locations_cbou)
only_in_tundra <- setdiff(unique_locations_tundra, unique_locations_cbou)
only_in_cbou <- setdiff(unique_locations_cbou, unique_locations_tundra)
# Output the results
cat("Common locations:", common_locations, "\n")## Common locations: BIO-TDN-018-01 BIO-TDN-018-02 BIO-TDN-018-03 BIO-TDN-018-04 BIO-TDN-018-05 BIO-TDN-019-01 BIO-TDN-019-02 BIO-TDN-019-03 BIO-TDN-019-04 BIO-TDN-019-05 BIO-TDN-028-01 BIO-TDN-028-02 BIO-TDN-028-03 BIO-TDN-028-05 BIO-TDN-028-06 BIO-TDN-031-01 BIO-TDN-031-02 BIO-TDN-031-03 BIO-TDN-031-04 BIO-TDN-031-05 BIO-TDN-034-01 BIO-TDN-034-02 BIO-TDN-034-03 BIO-TDN-034-04 BIO-TDN-034-05 BIO-TDN-034-NR BIO-TDN-036-01 BIO-TDN-036-02 BIO-TDN-036-03 BIO-TDN-036-04 BIO-TDN-036-05 BIO-TDN-037-01 BIO-TDN-037-03 BIO-TDN-037-04 BIO-TDN-037-05 BIO-TDN-037-06 BIO-TDN-042-01 BIO-TDN-042-02 BIO-TDN-042-03 BIO-TDN-042-06 BIO-TDN-042-07 BIO-TDN-044-01 BIO-TDN-044-02 BIO-TDN-044-03 BIO-TDN-044-07 BIO-TDN-044-08 BIO-TDN-045-01 BIO-TDN-045-02 BIO-TDN-045-03 BIO-TDN-045-05 BIO-TDN-045-06 BIO-TDN-046-01 BIO-TDN-046-02 BIO-TDN-046-03 BIO-TDN-046-04 BIO-TDN-046-05 BIO-TDN-046-NR BIO-TDN-049-01 BIO-TDN-049-02 BIO-TDN-049-03 BIO-TDN-049-04 BIO-TDN-049-05 BIO-TDN-051-01 BIO-TDN-051-02 BIO-TDN-051-03 BIO-TDN-051-05 BIO-TDN-051-06 BIO-TDN-052-02 BIO-TDN-052-04 BIO-TDN-052-05 BIO-TDN-052-06 BIO-TDN-052-07 BIO-TDN-053-01 BIO-TDN-053-02 BIO-TDN-053-03 BIO-TDN-053-04 BIO-TDN-053-05 BIO-TDN-054-01 BIO-TDN-054-02 BIO-TDN-054-04 BIO-TDN-054-05 BIO-TDN-054-08 BIO-TDN-055-01 BIO-TDN-055-02 BIO-TDN-055-04 BIO-TDN-055-05 BIO-TDN-055-06 BIO-TDN-056-01 BIO-TDN-056-02 BIO-TDN-056-03 BIO-TDN-056-04 BIO-TDN-056-05 BIO-TDN-057-01 BIO-TDN-057-02 BIO-TDN-057-03 BIO-TDN-057-05 BIO-TDN-057-06 BIO-TDN-058-01 BIO-TDN-058-02 BIO-TDN-058-03 BIO-TDN-058-05 BIO-TDN-058-06 BIO-TDN-059-01 BIO-TDN-059-02 BIO-TDN-059-03 BIO-TDN-059-04 BIO-TDN-059-05 BIO-TDN-060-01 BIO-TDN-060-02 BIO-TDN-060-04 BIO-TDN-060-05 BIO-TDN-060-06 BIO-TDN-061-01 BIO-TDN-061-02 BIO-TDN-061-03 BIO-TDN-061-05 BIO-TDN-061-06 BIO-TDN-062-01 BIO-TDN-062-02 BIO-TDN-062-03 BIO-TDN-062-04 BIO-TDN-062-05 BIO-TDN-164-01 BIO-TDN-164-02 BIO-TDN-164-04 BIO-TDN-164-05 BIO-TDN-164-06 BIO-TDN-165-01 BIO-TDN-165-02 BIO-TDN-165-04 BIO-TDN-165-05 BIO-TDN-165-06 BIO-TDN-165-NR BMS-CRU-004-01 BMS-CRU-004-02 BMS-CRU-004-03 BMS-CRU-004-04 BMS-CRU-004-05 BMS-CRU-007-02 BMS-CRU-007-03 BMS-CRU-007-05 BMS-CRU-007-06 BMS-CRU-007-08 BMS-CRU-008-01 BMS-CRU-008-02 BMS-CRU-008-03 BMS-CRU-008-04 BMS-CRU-008-05cat("Locations only in mod_dat_wk_tundra:", only_in_tundra, "\n")## Locations only in mod_dat_wk_tundra: BIO-TDN-009-01cat("Locations only in mod_dat_wk_cbou:", only_in_cbou, "\n")## Locations only in mod_dat_wk_cbou: BIO-TDN-050-01 BIO-TDN-050-03 BIO-TDN-050-04 BIO-TDN-050-06 BIO-TDN-050-07Map it out to see the cameras visually
# Create a leaflet map
map <- leaflet() %>%
addTiles() %>%
addCircleMarkers(data = mod_dat_wk_tundra[mod_dat_wk_tundra$location %in% common_locations,],
lat = ~latitude,
lng = ~longitude,
label = ~paste(location, " (Common)"),
color = "blue",
radius = 5) %>%
addCircleMarkers(data = mod_dat_wk_cbou[mod_dat_wk_cbou$location %in% common_locations,],
lat = ~latitude,
lng = ~longitude,
label = ~paste(location, " (Common)"),
color = "red",
radius = 5) %>%
addCircleMarkers(data = mod_dat_wk_tundra[mod_dat_wk_tundra$location %in% unique_locations_tundra,],
lat = ~latitude,
lng = ~longitude,
label = ~paste(location, " (Tundra)"),
color = "green",
radius = 5) %>%
addCircleMarkers(data = mod_dat_wk_cbou[mod_dat_wk_cbou$location %in% unique_locations_cbou,],
lat = ~latitude,
lng = ~longitude,
label = ~paste(location, " (Caribou)"),
color = "orange",
radius = 5) %>%
addLegend("bottomright",
colors = c("blue", "red", "green", "orange"),
labels = c("Common (Tundra)", "Common (Caribou)", "Unique (Tundra)", "Unique (Caribou)"))
# Print the map
mapVery annoyingly, there is one single camera with caribou detections which is not in a majority tundra-type landcover area.
BIO-TDN-050-06 has only two independent detections, both on 2021-11-24, first of a large bachelor group during the day and next of two females that night. Right in the transitional area between tundra and taiga but the camera micro site is very much needle leaf forest.
##Plotting
Local Scale Caribou
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk_cbou %>%
group_by(week_year, snow_season) %>%
summarise(total_caribou = sum(Caribou),total_effort = sum(n_days_effort))## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.#create capture rate
weekly_sum$total_weekly_cr<-((weekly_sum$total_caribou)/(weekly_sum$total_effort)*7)
# Plot the count data with color based on the 'season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = snow_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Caribou Count Summed Across Local Scale by Week",
x = "Week of the Year",
y = "Total Caribou Count") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8))## Warning: Removed 1 rows containing missing values (`geom_bar()`).
#plot the Capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = snow_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Caribou Weekly CR Across Local Scale by Week",
x = "Week of the Year",
y = "Total Caribou CR") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y.%m.%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8)) ## Warning: Removed 1 rows containing missing values (`geom_bar()`).
# Let's look at caribou detections summed across locations per week
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk_cbou %>%
group_by(week_year, four_season) %>%
summarise(total_caribou = sum(Caribou), total_effort = sum(n_days_effort))## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.# Create capture rate
weekly_sum$total_weekly_cr <- (weekly_sum$total_caribou / weekly_sum$total_effort) * 7
# Plot the count data with color based on the 'four_season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = four_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Caribou Count Summed Across Local Scale by Week",
x = "Week of the Year",
y = "Total Caribou Count") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8))## Warning: Removed 1 rows containing missing values (`geom_bar()`).
# Plot the capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = four_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Caribou Weekly CR Across Local Scale by Week",
x = "Week of the Year",
y = "Total Caribou CR") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) ## Warning: Removed 1 rows containing missing values (`geom_bar()`).
Muskox
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk_cbou %>%
group_by(week_year, snow_season) %>%
summarise(total_caribou = sum(Muskox),total_effort = sum(n_days_effort))## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.#create capture rate
weekly_sum$total_weekly_cr<-((weekly_sum$total_caribou)/(weekly_sum$total_effort)*7)
# Plot the count data with color based on the 'season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = snow_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Muskox Count Summed Across Local Scale by Week",
x = "Week of the Year",
y = "Total Muskox Count") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8))## Warning: Removed 1 rows containing missing values (`geom_bar()`).
#plot the Capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = snow_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Muskox Weekly CR Across Local Scale by Week",
x = "Week of the Year",
y = "Total Muskox CR") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y.%m.%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8)) ## Warning: Removed 1 rows containing missing values (`geom_bar()`).
# Let's look at caribou detections summed across locations per week
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk_cbou %>%
group_by(week_year, four_season) %>%
summarise(total_caribou = sum(Muskox), total_effort = sum(n_days_effort))## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.# Create capture rate
weekly_sum$total_weekly_cr <- (weekly_sum$total_caribou / weekly_sum$total_effort) * 7
# Plot the count data with color based on the 'four_season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = four_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Muskox Count Summed Across Local Scale by Week",
x = "Week of the Year",
y = "Total Muskox Count") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8))## Warning: Removed 1 rows containing missing values (`geom_bar()`).
# Plot the capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = four_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Muskox Weekly CR Across Local Scale by Week",
x = "Week of the Year",
y = "Total Muskox CR") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) ## Warning: Removed 1 rows containing missing values (`geom_bar()`).
Moose
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk_cbou %>%
group_by(week_year, snow_season) %>%
summarise(total_caribou = sum(Moose),total_effort = sum(n_days_effort))## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.#create capture rate
weekly_sum$total_weekly_cr<-((weekly_sum$total_caribou)/(weekly_sum$total_effort)*7)
# Plot the count data with color based on the 'season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = snow_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Moose Count Summed Across Local Scale by Week",
x = "Week of the Year",
y = "Total Moose Count") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8))## Warning: Removed 1 rows containing missing values (`geom_bar()`).
#plot the Capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = snow_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Moose Weekly CR Across Local Scale by Week",
x = "Week of the Year",
y = "Total Moose CR") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y.%m.%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8)) ## Warning: Removed 1 rows containing missing values (`geom_bar()`).
# Let's look at caribou detections summed across locations per week
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk_cbou %>%
group_by(week_year, four_season) %>%
summarise(total_caribou = sum(Moose), total_effort = sum(n_days_effort))## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.# Create capture rate
weekly_sum$total_weekly_cr <- (weekly_sum$total_caribou / weekly_sum$total_effort) * 7
# Plot the count data with color based on the 'four_season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = four_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Moose Count Summed Across Local Scale by Week",
x = "Week of the Year",
y = "Total Moose Count") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8))## Warning: Removed 1 rows containing missing values (`geom_bar()`).
# Plot the capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = four_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Moose Weekly CR Across Local Scale by Week",
x = "Week of the Year",
y = "Total Moose CR") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) ## Warning: Removed 1 rows containing missing values (`geom_bar()`).
Gray.Wolf
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk_cbou %>%
group_by(week_year, snow_season) %>%
summarise(total_caribou = sum(Gray.Wolf),total_effort = sum(n_days_effort))## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.#create capture rate
weekly_sum$total_weekly_cr<-((weekly_sum$total_caribou)/(weekly_sum$total_effort)*7)
# Plot the count data with color based on the 'season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = snow_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Gray Wolf Count Summed Across Local Scale by Week",
x = "Week of the Year",
y = "Total Gray Wolf Count") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8))## Warning: Removed 1 rows containing missing values (`geom_bar()`).
#plot the Capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = snow_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Gray Wolf Weekly CR Across Local Scale by Week",
x = "Week of the Year",
y = "Total Gray Wolf CR") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y.%m.%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8)) ## Warning: Removed 1 rows containing missing values (`geom_bar()`).
# Let's look at caribou detections summed across locations per week
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk_cbou %>%
group_by(week_year, four_season) %>%
summarise(total_caribou = sum(Gray.Wolf), total_effort = sum(n_days_effort))## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.# Create capture rate
weekly_sum$total_weekly_cr <- (weekly_sum$total_caribou / weekly_sum$total_effort) * 7
# Plot the count data with color based on the 'four_season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = four_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Gray Wolf Count Summed Across Local Scale by Week",
x = "Week of the Year",
y = "Total Gray Wolf Count") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8))## Warning: Removed 1 rows containing missing values (`geom_bar()`).
# Plot the capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = four_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Gray Wolf Weekly CR Across Local Scale by Week",
x = "Week of the Year",
y = "Total Gray Wolf CR") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) ## Warning: Removed 1 rows containing missing values (`geom_bar()`).
Grizzly.Bear
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk_cbou %>%
group_by(week_year, snow_season) %>%
summarise(total_caribou = sum(Grizzly.Bear),total_effort = sum(n_days_effort))## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.#create capture rate
weekly_sum$total_weekly_cr<-((weekly_sum$total_caribou)/(weekly_sum$total_effort)*7)
# Plot the count data with color based on the 'season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = snow_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Grizzly Bear Count Summed Across Local Scale by Week",
x = "Week of the Year",
y = "Total Grizzly Bear Count") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8))## Warning: Removed 1 rows containing missing values (`geom_bar()`).
#plot the Capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = snow_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Grizzly Bear Weekly CR Across Local Scale by Week",
x = "Week of the Year",
y = "Total Grizzly Bear CR") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y.%m.%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size=8)) ## Warning: Removed 1 rows containing missing values (`geom_bar()`).
# Let's look at caribou detections summed across locations per week
# Group by week_year and summarize the total Caribou count and effort
weekly_sum <- mod_dat_wk_cbou %>%
group_by(week_year, four_season) %>%
summarise(total_caribou = sum(Grizzly.Bear), total_effort = sum(n_days_effort))## `summarise()` has grouped output by 'week_year'. You can override using the
## `.groups` argument.# Create capture rate
weekly_sum$total_weekly_cr <- (weekly_sum$total_caribou / weekly_sum$total_effort) * 7
# Plot the count data with color based on the 'four_season' column
ggplot(weekly_sum, aes(x = week_year, y = total_caribou, fill = four_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Grizzly Bear Count Summed Across Local Scale by Week",
x = "Week of the Year",
y = "Total Grizzly Bear Count") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8))## Warning: Removed 1 rows containing missing values (`geom_bar()`).
# Plot the capture rate
ggplot(weekly_sum, aes(x = week_year, y = total_weekly_cr, fill = four_season)) +
geom_bar(stat = "identity", position = "dodge", color = "black") +
labs(title = "Total Grizzly Bear Weekly CR Across Local Scale by Week",
x = "Week of the Year",
y = "Total Grizzly Bear CR") +
theme_minimal() +
scale_x_date(labels = scales::date_format("%Y-%m-%d"), date_breaks = "1 week") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 8)) ## Warning: Removed 1 rows containing missing values (`geom_bar()`).
Local, Multiseason, WEEK
##Target analysis
# Reorder the levels of the target_cons column
mod_dat_wk_cbou$target_cons <- fct_relevel(mod_dat_wk_cbou$target_cons, "None")
#check the levels were reordered correctly. None should be first to act as the control feature.
levels(mod_dat_wk_cbou$target_cons)## [1] "None" "Esker" "Shoreline" "Trail" "Wetland"# Group by target_cons and count unique locations within each group
unique_cbou_locations_by_target <- mod_dat_wk_cbou %>%
group_by(target_cons) %>%
summarise(unique_locations = n_distinct(location))
# Print the result
print(unique_cbou_locations_by_target)## # A tibble: 5 × 2
## target_cons unique_locations
## <fct> <int>
## 1 None 108
## 2 Esker 7
## 3 Shoreline 13
## 4 Trail 5
## 5 Wetland 20#run anova to determine whether targets are actually significantly related to our response variable (Caribou occurrences in a given week at a given camera)
region_wk_cbou_aov <- aov(Caribou ~ target_cons, data=mod_dat_wk_cbou)
summary(region_wk_cbou_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_cons 4 26 6.576 2.479 0.042 *
## Residuals 6351 16844 2.652
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_cbou_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_cons, data = mod_dat_wk_cbou)
##
## $target_cons
## diff lwr upr p adj
## Esker-None 0.14286874 -0.2282522 0.51398970 0.8317097
## Shoreline-None 0.02273957 -0.1869334 0.23241255 0.9983316
## Trail-None -0.19579310 -0.4866200 0.09503386 0.3522669
## Wetland-None -0.13667776 -0.2998466 0.02649109 0.1496063
## Shoreline-Esker -0.12012917 -0.5361447 0.29588642 0.9342055
## Trail-Esker -0.33866183 -0.8009242 0.12360057 0.2665033
## Wetland-Esker -0.27954650 -0.6741733 0.11508027 0.2998968
## Trail-Shoreline -0.21853267 -0.5648252 0.12775985 0.4204316
## Wetland-Shoreline -0.15941733 -0.4083395 0.08950488 0.4048496
## Wetland-Trail 0.05911534 -0.2611655 0.37939619 0.9870349#run anova on hydro targets
region_wk_cbou_hydro_aov <- aov(Caribou ~ target_hydro, data=mod_dat_wk_cbou)
summary(region_wk_cbou_hydro_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_hydro 2 14 6.999 2.638 0.0716 .
## Residuals 6353 16857 2.653
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_cbou_hydro_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_hydro, data = mod_dat_wk_cbou)
##
## $target_hydro
## diff lwr upr p adj
## Shoreline-None 0.02817302 -0.1513042 0.207650247 0.9280829
## Wetland-None -0.13124431 -0.2705580 0.008069353 0.0698140
## Wetland-Shoreline -0.15941733 -0.3733242 0.054489516 0.1878788Not significant (wetland-none is close though at 0.07)
#run anova on linear targets
region_wk_cbou_linear_aov <- aov(Caribou ~ target_linear, data=mod_dat_wk_cbou)
summary(region_wk_cbou_linear_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_linear 2 12 5.765 2.172 0.114
## Residuals 6353 16859 2.654# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_cbou_linear_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_linear, data = mod_dat_wk_cbou)
##
## $target_linear
## diff lwr upr p adj
## Esker-None 0.1612578 -0.1565214 0.47903699 0.4594069
## Trail-None -0.1774040 -0.4258569 0.07104883 0.2153357
## Trail-Esker -0.3386618 -0.7359278 0.05860418 0.1126497Not Significant
Local, Summer, WEEKLY detection data
# Create a new data frame for summer
mod_dat_wk_cbou_summer <- mod_dat_wk_cbou[mod_dat_wk_cbou$four_season == "Summer", ]##Target analysis
# Group by target_cons and count unique locations within each group
unique_cbou_locations_by_target_smr <- mod_dat_wk_cbou_summer %>%
group_by(target_cons) %>%
summarise(unique_locations = n_distinct(location))
# Print the result
print(unique_cbou_locations_by_target_smr)## # A tibble: 5 × 2
## target_cons unique_locations
## <fct> <int>
## 1 None 108
## 2 Esker 7
## 3 Shoreline 13
## 4 Trail 5
## 5 Wetland 20#run anova to see if target is important
region_wk_cbou_smr_aov <- aov(Caribou ~ target_cons, data=mod_dat_wk_cbou_summer)
summary(region_wk_cbou_smr_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_cons 4 14 3.479 2.453 0.0441 *
## Residuals 2272 3223 1.419
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Perform Tukey's HSD test to interpret ANOVA
tukey_results <- TukeyHSD(region_wk_cbou_smr_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_cons, data = mod_dat_wk_cbou_summer)
##
## $target_cons
## diff lwr upr p adj
## Esker-None -0.021674699 -0.47140773 0.42805833 0.9999326
## Shoreline-None 0.188071890 -0.06114205 0.43728583 0.2379321
## Trail-None -0.025378402 -0.37742875 0.32667194 0.9996660
## Wetland-None -0.153482977 -0.35599721 0.04903125 0.2339582
## Shoreline-Esker 0.209746589 -0.29170087 0.71119404 0.7841851
## Trail-Esker -0.003703704 -0.56341917 0.55601177 1.0000000
## Wetland-Esker -0.131808279 -0.61176033 0.34814378 0.9446105
## Trail-Shoreline -0.213450292 -0.62953826 0.20263768 0.6274293
## Wetland-Shoreline -0.341554868 -0.64189096 -0.04121878 0.0165115
## Wetland-Trail -0.128104575 -0.51801947 0.26181032 0.8981167#run anova on hydro targets
region_wk_cbou_smr_hydro_aov <- aov(Caribou ~ target_hydro, data=mod_dat_wk_cbou_summer)
summary(region_wk_cbou_smr_hydro_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_hydro 2 14 6.920 4.882 0.00766 **
## Residuals 2274 3223 1.417
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_cbou_smr_hydro_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_hydro, data = mod_dat_wk_cbou_summer)
##
## $target_hydro
## diff lwr upr p adj
## Shoreline-None 0.1900115 -0.02308802 0.40311107 0.0918697
## Wetland-None -0.1515433 -0.32433252 0.02124583 0.0992151
## Wetland-Shoreline -0.3415549 -0.59945484 -0.08365489 0.0054503 Df Sum Sq Mean Sq F value Pr(>F) target_hydro 2 14 6.912 4.952 0.00715 **
diff lwr upr p adjWetland-Shoreline -0.3415549 -0.59747378 -0.08563595 0.0050335
******Not significantly different from NONE
#run anova on linear targets
region_wk_cbou_smr_linear_aov <- aov(Caribou ~ target_linear, data=mod_dat_wk_cbou_summer)
summary(region_wk_cbou_smr_linear_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_linear 2 0 0.0239 0.017 0.983
## Residuals 2274 3237 1.4236# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_cbou_smr_linear_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_linear, data = mod_dat_wk_cbou_summer)
##
## $target_linear
## diff lwr upr p adj
## Esker-None -0.016408814 -0.4019824 0.3691648 0.9945233
## Trail-None -0.020112518 -0.3212246 0.2809996 0.9865630
## Trail-Esker -0.003703704 -0.4853612 0.4779538 0.9998207Not Significant
Local, Fall, WEEKLY detection data
# Create a new data frame for winter
mod_dat_wk_cbou_fall <- mod_dat_wk_cbou[mod_dat_wk_cbou$four_season == "Fall", ]##Target analysis
# Group by target_cons and count unique locations within each group
unique_cbou_locations_by_target_fall <- mod_dat_wk_cbou_fall %>%
group_by(target_cons) %>%
summarise(unique_locations = n_distinct(location))
# Print the result
print(unique_cbou_locations_by_target_fall)## # A tibble: 5 × 2
## target_cons unique_locations
## <fct> <int>
## 1 None 104
## 2 Esker 4
## 3 Shoreline 13
## 4 Trail 5
## 5 Wetland 20#run anova to see if target is important
region_wk_cbou_fll_aov <- aov(Caribou ~ target_cons, data=mod_dat_wk_cbou_fall)
summary(region_wk_cbou_fll_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_cons 4 58 14.49 1.403 0.231
## Residuals 1014 10470 10.33# Perform Tukey's HSD test to interpret ANOVA
tukey_results <- TukeyHSD(region_wk_cbou_fll_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_cons, data = mod_dat_wk_cbou_fall)
##
## $target_cons
## diff lwr upr p adj
## Esker-None 0.1494062 -1.5712267 1.8700391 0.9992990
## Shoreline-None -0.1383431 -1.1839603 0.9072742 0.9963502
## Trail-None -1.0870610 -2.5299283 0.3558062 0.2390066
## Wetland-None -0.3787046 -1.1907952 0.4333859 0.7071021
## Shoreline-Esker -0.2877493 -2.2484541 1.6729555 0.9945514
## Trail-Esker -1.2364672 -3.4348561 0.9619217 0.5385550
## Wetland-Esker -0.5281108 -2.3748753 1.3186536 0.9359919
## Trail-Shoreline -0.9487179 -2.6708247 0.7733888 0.5591328
## Wetland-Shoreline -0.2403616 -1.4826485 1.0019253 0.9844157
## Wetland-Trail 0.7083564 -0.8828153 2.2995281 0.7417760#run anova on hydro targets
region_wk_cbou_fll_hydro_aov <- aov(Caribou ~ target_hydro, data=mod_dat_wk_cbou_fall)
summary(region_wk_cbou_fll_hydro_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_hydro 2 13 6.547 0.633 0.531
## Residuals 1016 10515 10.349# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_cbou_fll_hydro_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_hydro, data = mod_dat_wk_cbou_fall)
##
## $target_hydro
## diff lwr upr p adj
## Shoreline-None -0.0905109 -0.9860893 0.8050675 0.9694576
## Wetland-None -0.3308725 -1.0246129 0.3628680 0.5022606
## Wetland-Shoreline -0.2403616 -1.3086092 0.8278861 0.8575650Not Significant
#run anova on linear targets
region_wk_cbou_fll_linear_aov <- aov(Caribou ~ target_linear, data=mod_dat_wk_cbou_fall)
summary(region_wk_cbou_fll_linear_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_linear 2 41 20.37 1.974 0.139
## Residuals 1016 10487 10.32# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_cbou_fll_linear_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_linear, data = mod_dat_wk_cbou_fall)
##
## $target_linear
## diff lwr upr p adj
## Esker-None 0.2159652 -1.255703 1.6876329 0.9367034
## Trail-None -1.0205021 -2.252478 0.2114736 0.1269366
## Trail-Esker -1.2364672 -3.124384 0.6514495 0.2739310Not Significant
Local, Winter, WEEKLY detection data
# Create a new data frame for winter
mod_dat_wk_cbou_winter <- mod_dat_wk_cbou[mod_dat_wk_cbou$four_season == "Winter", ]##Target analysis
# Group by target_cons and count unique locations within each group
unique_cbou_locations_by_target_wtr <- mod_dat_wk_cbou_winter %>%
group_by(target_cons) %>%
summarise(unique_locations = n_distinct(location))
# Print the result
print(unique_cbou_locations_by_target_wtr)## # A tibble: 5 × 2
## target_cons unique_locations
## <fct> <int>
## 1 None 101
## 2 Esker 3
## 3 Shoreline 13
## 4 Trail 5
## 5 Wetland 19#run anova to see if target is important
region_wk_cbou_wtr_aov <- aov(Caribou ~ target_cons, data=mod_dat_wk_cbou_winter)
summary(region_wk_cbou_wtr_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_cons 4 7.3 1.8336 4.833 0.00071 ***
## Residuals 1551 588.4 0.3794
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Perform Tukey's HSD test to interpret ANOVA
tukey_results <- TukeyHSD(region_wk_cbou_wtr_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_cons, data = mod_dat_wk_cbou_winter)
##
## $target_cons
## diff lwr upr p adj
## Esker-None 0.378046595 0.1073555 0.64873769 0.0013369
## Shoreline-None -0.117827191 -0.2910513 0.05539690 0.3409280
## Trail-None -0.017921147 -0.2374046 0.20156234 0.9994520
## Wetland-None -0.027804469 -0.1485355 0.09292653 0.9704193
## Shoreline-Esker -0.495873786 -0.8092576 -0.18248998 0.0001603
## Trail-Esker -0.395967742 -0.7371068 -0.05482873 0.0134634
## Wetland-Esker -0.405851064 -0.6935641 -0.11813801 0.0011496
## Trail-Shoreline 0.099906044 -0.1704802 0.37029233 0.8512838
## Wetland-Shoreline 0.090022723 -0.1087528 0.28879827 0.7297651
## Wetland-Trail -0.009883322 -0.2500460 0.23027940 0.9999640#run anova on hydro targets
region_wk_cbou_wtr_hydro_aov <- aov(Caribou ~ target_hydro, data=mod_dat_wk_cbou_winter)
summary(region_wk_cbou_wtr_hydro_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_hydro 2 1.8 0.8794 2.299 0.101
## Residuals 1553 594.0 0.3825# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_cbou_wtr_hydro_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_hydro, data = mod_dat_wk_cbou_winter)
##
## $target_hydro
## diff lwr upr p adj
## Shoreline-None -0.12933027 -0.27820950 0.01954896 0.1035810
## Wetland-None -0.03930755 -0.14267890 0.06406380 0.6454035
## Wetland-Shoreline 0.09002272 -0.08142462 0.26147007 0.4345380Not Significant
#run anova on linear targets
region_wk_cbou_wtr_linear_aov <- aov(Caribou ~ target_linear, data=mod_dat_wk_cbou_winter)
summary(region_wk_cbou_wtr_linear_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_linear 2 6.0 2.9812 7.85 0.000405 ***
## Residuals 1553 589.8 0.3798
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_cbou_wtr_linear_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_linear, data = mod_dat_wk_cbou_winter)
##
## $target_linear
## diff lwr upr p adj
## Esker-None 0.390887208 0.1591766 0.6225978 0.0002334
## Trail-None -0.005080534 -0.1925603 0.1823992 0.9977743
## Trail-Esker -0.395967742 -0.6891631 -0.1027724 0.0044555 Df Sum Sq Mean Sq F value Pr(>F) target_linear 2 6.0 2.9911 7.94 0.000371 ***
diff lwr upr p adjEsker-None 0.391574966 0.1607959 0.6223541 0.0002123
Local, Spring, WEEKLY detection data
# Create a new data frame for spring
mod_dat_wk_cbou_spri <- mod_dat_wk_cbou[mod_dat_wk_cbou$four_season == "Spring", ]##Target analysis
# Group by target_cons and count unique locations within each group
unique_cbou_locations_by_target_spr <- mod_dat_wk_cbou_spri %>%
group_by(target_cons) %>%
summarise(unique_locations = n_distinct(location))
# Print the result
print(unique_cbou_locations_by_target_spr)## # A tibble: 5 × 2
## target_cons unique_locations
## <fct> <int>
## 1 None 102
## 2 Esker 3
## 3 Shoreline 13
## 4 Trail 5
## 5 Wetland 19#run anova to see if target is important
region_wk_cbou_spr_aov <- aov(Caribou ~ target_cons, data=mod_dat_wk_cbou_spri)
summary(region_wk_cbou_spr_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_cons 4 0.7 0.1677 0.383 0.821
## Residuals 1499 656.9 0.4382# Perform Tukey's HSD test to interpret ANOVA
tukey_results <- TukeyHSD(region_wk_cbou_spr_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_cons, data = mod_dat_wk_cbou_spri)
##
## $target_cons
## diff lwr upr p adj
## Esker-None -0.08715596 -0.4393827 0.26507074 0.9616486
## Shoreline-None 0.01520624 -0.1543138 0.18472626 0.9992054
## Trail-None -0.03261051 -0.2824718 0.21725078 0.9965482
## Wetland-None -0.04813157 -0.1857687 0.08950553 0.8750307
## Shoreline-Esker 0.10236220 -0.2807866 0.48551098 0.9496320
## Trail-Esker 0.05454545 -0.3703035 0.47939438 0.9967615
## Wetland-Esker 0.03902439 -0.3311241 0.40917290 0.9984998
## Trail-Shoreline -0.04781675 -0.3396557 0.24402224 0.9917116
## Wetland-Shoreline -0.06333781 -0.2675027 0.14082708 0.9156957
## Wetland-Trail -0.01552106 -0.2900697 0.25902753 0.9998725#run anova on hydro targets
region_wk_cbou_spr_hydro_aov <- aov(Caribou ~ target_hydro, data=mod_dat_wk_cbou_spri)
summary(region_wk_cbou_spr_hydro_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_hydro 2 0.4 0.2109 0.482 0.618
## Residuals 1501 657.2 0.4378# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_cbou_spr_hydro_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_hydro, data = mod_dat_wk_cbou_spri)
##
## $target_hydro
## diff lwr upr p adj
## Shoreline-None 0.01874446 -0.1262726 0.16376147 0.9505772
## Wetland-None -0.04459336 -0.1621119 0.07292516 0.6465706
## Wetland-Shoreline -0.06333781 -0.2386331 0.11195750 0.6733747Not Significant
#run anova on linear targets
region_wk_cbou_spr_linear_aov <- aov(Caribou ~ target_linear, data=mod_dat_wk_cbou_spri)
summary(region_wk_cbou_spr_linear_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## target_linear 2 0.2 0.1054 0.241 0.786
## Residuals 1501 657.4 0.4380# Perform Tukey's HSD test
tukey_results <- TukeyHSD(region_wk_cbou_spr_linear_aov)
# View the results
print(tukey_results)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Caribou ~ target_linear, data = mod_dat_wk_cbou_spri)
##
## $target_linear
## diff lwr upr p adj
## Esker-None -0.08157525 -0.3831898 0.2200393 0.8011177
## Trail-None -0.02702979 -0.2403878 0.1863282 0.9524751
## Trail-Esker 0.05454545 -0.3102870 0.4193779 0.9344405Not Significant
#SUMMARIZE
Create tables summarizing effort and independent detections for each focal species in each subset dataset.
# List of dataset names
dataset_names <- c("mod_dat_wk", "mod_dat_wk_summer", "mod_dat_wk_fall", "mod_dat_wk_winter", "mod_dat_wk_spri",
"mod_dat_wk_cbou", "mod_dat_wk_cbou_summer", "mod_dat_wk_cbou_fall", "mod_dat_wk_cbou_winter",
"mod_dat_wk_cbou_spri")
# Get the dataframes from the environment
datasets <- mget(dataset_names)
# Function to create summary dataframe for each dataset
summarize_dataset <- function(df, dataset_name) {
summary_df <- data.frame(
Dataset = dataset_name,
Total_Effort = sum(df$n_days_effort),
Caribou_Detections = sum(df$Caribou, na.rm = TRUE),
Muskox_Detections = sum(df$Muskox, na.rm = TRUE),
Moose_Detections = sum(df$Moose, na.rm = TRUE),
Gray_Wolf_Detections = sum(df$Gray.Wolf, na.rm = TRUE),
Grizzly_Bear_Detections = sum(df$Grizzly.Bear, na.rm = TRUE)
)
return(summary_df)
}
# Create summary dataframes for each dataset
summary_dfs <- mapply(summarize_dataset, datasets, dataset_names, SIMPLIFY = FALSE)
# Combine summary dataframes into a single dataframe
summary_all <- do.call(rbind, summary_dfs)
# Print summary dataframe
print(summary_all)## Dataset Total_Effort Caribou_Detections
## mod_dat_wk mod_dat_wk 83309 2348
## mod_dat_wk_summer mod_dat_wk_summer 30002 372
## mod_dat_wk_fall mod_dat_wk_fall 12490 1633
## mod_dat_wk_winter mod_dat_wk_winter 19238 224
## mod_dat_wk_spri mod_dat_wk_spri 21579 119
## mod_dat_wk_cbou mod_dat_wk_cbou 40039 2348
## mod_dat_wk_cbou_summer mod_dat_wk_cbou_summer 14733 372
## mod_dat_wk_cbou_fall mod_dat_wk_cbou_fall 5971 1633
## mod_dat_wk_cbou_winter mod_dat_wk_cbou_winter 9226 224
## mod_dat_wk_cbou_spri mod_dat_wk_cbou_spri 10109 119
## Muskox_Detections Moose_Detections Gray_Wolf_Detections
## mod_dat_wk 746 180 130
## mod_dat_wk_summer 599 97 42
## mod_dat_wk_fall 54 13 34
## mod_dat_wk_winter 22 22 15
## mod_dat_wk_spri 71 48 39
## mod_dat_wk_cbou 212 29 57
## mod_dat_wk_cbou_summer 180 19 6
## mod_dat_wk_cbou_fall 11 2 27
## mod_dat_wk_cbou_winter 3 1 13
## mod_dat_wk_cbou_spri 18 7 11
## Grizzly_Bear_Detections
## mod_dat_wk 38
## mod_dat_wk_summer 25
## mod_dat_wk_fall 0
## mod_dat_wk_winter 0
## mod_dat_wk_spri 13
## mod_dat_wk_cbou 37
## mod_dat_wk_cbou_summer 25
## mod_dat_wk_cbou_fall 0
## mod_dat_wk_cbou_winter 0
## mod_dat_wk_cbou_spri 12#save table
write.csv(summary_all, file = "Tables/GLMM_DataSet_CountandEffort_Comparison.csv", row.names = F) Dataset Total_Effort Caribou_Detections Muskox_Detections Moose_Detectionsmod_dat_wk mod_dat_wk 83309 2348 746 180 mod_dat_wk_summer mod_dat_wk_summer 30002 372 599 97 mod_dat_wk_fall mod_dat_wk_fall 12490 1633 54 13 mod_dat_wk_winter mod_dat_wk_winter 19238 224 22 22 mod_dat_wk_spri mod_dat_wk_spri 21579 119 71 48 mod_dat_wk_cbou mod_dat_wk_cbou 40039 2348 212 29 mod_dat_wk_cbou_summer mod_dat_wk_cbou_summer 14733 372 180 19 mod_dat_wk_cbou_fall mod_dat_wk_cbou_fall 5971 1633 11 2 mod_dat_wk_cbou_winter mod_dat_wk_cbou_winter 9226 224 3 1 mod_dat_wk_cbou_spri mod_dat_wk_cbou_spri 10109 119 18 7 Gray_Wolf_Detections Grizzly_Bear_Detections mod_dat_wk 130 38 mod_dat_wk_summer 42 25 mod_dat_wk_fall 34 0 mod_dat_wk_winter 15 0 mod_dat_wk_spri 39 13 mod_dat_wk_cbou 57 37 mod_dat_wk_cbou_summer 6 25 mod_dat_wk_cbou_fall 27 0 mod_dat_wk_cbou_winter 13 0 mod_dat_wk_cbou_spri 11 12
#SAVE DATA AS CSV
##Regional data sets
#multiseason
#write.csv(mod_dat_wk, file="Processed_Data/glmm_wk_region_multi.csv", row.names = F)
#summer
#write.csv(mod_dat_wk_summer, file="Processed_Data/glmm_wk_region_smr.csv", row.names = F)
#fall
#write.csv(mod_dat_wk_fall, file="Processed_Data/glmm_wk_region_fall.csv", row.names = F)
#winter
#write.csv(mod_dat_wk_winter, file="Processed_Data/glmm_wk_region_wtr.csv", row.names = F)
#spring
#write.csv(mod_dat_wk_spri, file="Processed_Data/glmm_wk_region_spr.csv", row.names = F)##Local data sets
#multiseason
#write.csv(mod_dat_wk_cbou, file="Processed_Data/glmm_wk_local_multi.csv", row.names = F)
#summer
#write.csv(mod_dat_wk_cbou_summer, file="Processed_Data/glmm_wk_local_smr.csv", row.names = F)
#fall
#write.csv(mod_dat_wk_cbou_fall, file="Processed_Data/glmm_wk_local_fall.csv", row.names = F)
#winter
#write.csv(mod_dat_wk_cbou_winter, file="Processed_Data/glmm_wk_local_wtr.csv", row.names = F)
#spring
#write.csv(mod_dat_wk_cbou_spri, file="Processed_Data/glmm_wk_local_spr.csv", row.names = F)#regional only dataset
#mod_dat_wk_region_only <- subset(mod_dat_wk, spatial_scale == "regional")